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ABSTRACT 


A  binary  collision  simulation  of  the  radiation  damage  problem  is 
proposed.  Energy  thresholds  and  parameter  adjustment  are  discussed. 
Preliminary  results  show  strong  focusing  in  the  (110)  direction  and 
weaker  focusing  in  the  (100)  direction. 

Although  not  complete,  the  model  is  expected,  with  suggested 
modifications,  to  adequately  represent  the  radiation  damage  problem. 
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1 .  Introduction 


The  interaction  of  energetic  particles  with  matter  may  be  divided 
into  two  energetic  regions  of  interest: 

1)  High  energy  events  in  which  interactions  are  primarily 
inelastic  electronic  processes  of  excitation  and  ionization, 
and 

2)  Low  energy  events  in  which  the  interactions  are  essentially 
elastic  interactions  between  whole  atoms. 

There  is  an  intermediate  energy  range  in  which  both  processes  are  im¬ 
portant,  but  until  the  two  limiting  cases  are  better  understood,  attempts 

1* 

to  understand  the  intermediate  range  seem  futile.  The  threshold  energy 
for  which  electronic  excitation  and  ionization  processes  become  impor¬ 
tant  depends  upon  the  mass  of  the  energetic  particle.  An  approximate 
formula  for  this  threshold  energy  is  E  =  A  kev.  ,  where  A  is  the  atomic 
mass  of  the  energetic  particle.  The  threshold  for  electronic  processes 

seems  to  be  relatively  insensitive  to  the  atomic  mass  of  the  target 
2 

material.  For  metals  this  threshold  is  of  the  order  of  10-100  kev.  The 
high  energy  processes  of  electronic  excitation  and  ionization  have  re¬ 
ceived  considerable  attention  and  are  reasonably  well  understood. 

We  will  confine  our  attention  to  the  energy  region  below  this  threshold 
for  inelastic  electronic  processes. 

Many  physical  properties  of  a  material  are  changed  when  the 
7 

material  is  irradiated.  The  Young's  modulus  and  the  electrical 
*  All  footnotes  refer  to  the  bibliography. 
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resistivity  are  among  those  properties  of  materials  which  show  most 
significant  changes.,  For  metals  and  semiconductors,  the  residual 
electrical  resistivity  is  most  sensitive  to  radiation. 

Semiconductors,  for  example,  may  show  changes  in  electrical 

7 

resistivity  of  up  to  four  orders  of  magnitude,  but  there  is  no  known 
simple  relation  between  the  magnitude  of  the  change  in  electrical 
resistivity  and  the  amount  of  radiative  energy  absorbed  by  a  semi¬ 
conductor.  Radiation  changes  some  semiconductors  from  n-type  to  the 
p-type . 

The  radiation  damage  problem  involves  chemical  bonding  between 
the  atoms  of  the  target  material  as  well  as  elastic  collisions  between 
atoms.  For  low  energy  collisions,  chemical  bonding  forces  vary  widely 
for  different  atomic  species,  but  are  of  the  same  order  of  magnitude  as 
the  elastic  repulsive  forces.  Therefore,  to  reduce  the  complexity  of 
the  problem,  we  have  limited  our  investigation  to  simple  crystals  of 
face  centered  cubic  metals,  with  special  emphasis  on  copper. 

Davies,  et  al.  have  developed  a  chemical  technique  for  removing 

g 

thin  layers  (2>  37  A)  from  aluminum  foils.  The  thickness  of  the  layers 
is  easily  controlled.  The  technique  has  been  widely  used  to  study  the 
range  of  radioactive  ions  in  aluminum.^  Ranges  have  been  measured 
for  Cs137,  Rb^S,  Na^,  and  K^.  For  example,  the  mean  range  of 
10.5  kev.  Na^4  ions  in  aluminum  was  found  to  be  about  250  Angstroms. 

When  a  metal  is  irradiated  with  heavy  energetic  ions  ( >  50  kev.)  a 
significant  number  of  atoms  are  ejected  from  the  target  through  the  surface 
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which  is  bombarded  by  the  incident  ion  beam.  This  phenomenon  is 
known  as  sputtering.  Under  suitable  conditions  as  many  as  40  atoms  or 
ions  may  be  sputtered  for  each  incident  particle.  ^  Examination  of  the 
angular  distribution  of  sputtered  particles  reveals  definite  preferential 

ejection  in  directions  which  correspond  to  (110),  (100),  and  (ill)  lattice 

,  12-16 
directions . 
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2.  Lattice  Effects 


Whether  a  lattice  atom  leaves  its  original  site  or  not  after  a 
collision  depends  upon  the  energy  absorbed  by  the  atom  and  its  final 
direction  of  motion  with  respect  to  the  lattice.  If  an  energetic  lattice 
atom  does  not  leave  its  site,  it  will  transfer  energy  to  its  neighbors 
through  elastic  collisions  until  thermal  equilibrium  is  reached.  Be¬ 
cause  of  crystalline  structure,  the  dissipation  of  energy  is  expected  to 
be  anisotropic.  The  only  exception  to  the  eventual  thermalization  of 
the  excess  energy  of  a  bound  lattice  atom  is  the  case  in  which  sufficient 
energy  reaches  a  surface  of  the  material  to  knock  one  or  more  atoms  away 
from  the  surface . 

If  a  lattice  atom  is  given  sufficient  energy  to  leave  its  site  per¬ 
manently,  and  no  other  atom  replaces  the  energetic  atom,  a  vacant 
lattice  site  is  created.  The  surrounding  lattice  atoms  relax  slightly 

toward  the  vacant  site,  but  the  lattice  structure  remains  essentially  the 

17 

same  as  that  of  a  perfect  lattice.  The  energetic  atom  travels  through 

the  lattice  giving  up  energy  through  elastic  collisions  with  the  atoms 

near  its  path  until  it  no  longer  has  sufficient  energy  to  force  its  way 

between  the  atoms  of  the  crystal.  If  there  are  no  vacant  lattice  sites 

near  the  position  at  which  an  energetic  atom  is  brought  to  rest,  it  will 

be  forced  to  occupy  some  position  other  than  a  regular  lattice  site. 

Three  possible  equilibrium  positions  have  been  suggested  for  these 

interstitial  atoms  in  a  face  centered  cubic  crystal:  the  crowdion ,  the 

17 

split  interstitial  pair,  and  the  center  of  a  unit  cell.  The  crowdion 
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consists  of  an  extra  atom  in  a  (110)  line  of  atoms.  In  the  split  inter¬ 
stitial  pair  configuration,  the  interstitial  shares  a  regular  lattice  site 
with  another  lattice  atom  so  that  each  is  symmetrically  displaced  from 
the  equilibrium  position  in  a  (100)  direction.  An  atom  in  the  center  of 
a  unit  cell  may  also  be  described  as  a  crowdion  in  a  (100)  direction. 

The  lattice  atoms  surrounding  an  interstitial  are  displaced  from  their 
regular  lattice  sites  to  accommodate  the  extra  atom.  More  will  be  said 
about  the  interstitial  configuration  later. 

An  energetic  atom  (knock-on)  may  transfer  essentially  all  of  its 

18 

energy  to  a  single  lattice  atom  in  a  head-on  collision.  In  this  case, 
the  lattice  atom  becomes  the  knock-on,  leaving  a  site  for  the  former 
knock-on  to  occupy.  The  new  knock-on  may,  in  turn,  transfer  essentially 
all  of  its  energy  to  another  atom.  Thus,  a  chain  of  replacements  may 
transport  an  interstitial  to  a  position  which  is  several  times  the  dis¬ 
tance  between  nearest  neighbors  from  the  corresponding  vacancy.  Re¬ 
placement  chains  may  be  cyclic  so  that  the  last  replacement  in  the  chain 

17 

fills  the  initial  vacancy,  leaving  no  vacancies  or  interstitials.  The 
configuration  consisting  of  a  vacancy  and  a  nearby  interstitial,  known 
as  a  Frenkel  pair,  may  be  stable  or  unstable.  The  stability  of  a  Frenkel 
pair  depends  upon  both  the  distance  of  separation  and  the  orientation 

(with  respect  to  the  lattice)  of  the  interstitial  and  the  vacancy. 

19 

Silsbee  has  shown  that  an  isolated  line  of  hard  spheres  may 
focus  and  transfer  energy  along  the  line.  The  condition  for  focusing 
is  that  s/d<.  cos  ^  ,  where  s  is  the  distance  between  adjacent  spheres, 
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d  is  the  diameter  of  the  spheres ,  and  @  is  the  angle  between  the 
velocity  of  the  first  atom  in  the  line  and  the  line.  The  hard  core 
approximation  (diameter  of  spheres  equal  to  distance  of  closest 
approach  for  a  head-on  collision)  is  useful  in  this  work.  For  reason¬ 
able  repulsive  potentials  such  as  the  Coulomb,  Bohr,  Born-Mayer,  etc.  , 
this  modified  model  predicts  a  focusing  effect  below  a  threshold  energy, 
and  a  de-focusing  effect  above  this  threshold  energy.^  When  the 
Gibson  potential  #2  simulates  copper,  this  model  predicts  focusing  in 
the  (110)  directions  below  about  60  ev.  ,  with  sharp  focusing  below 
about  15  ev.  The  presence  of  adjacent  lines  of  atoms  will  enhance 
the  focusing  effect  found  in  an  isolated  line  of  atoms .  Although  the 
hard  sphere  models  used  to  demonstrate  the  phenomenon  of  energy 
focusing  contain  rather  drastic  assumptions,  the  same  sort  of  effect  is 

expected  from  models  based  on  more  realistic  potentials. 

20 

A  channeling  effect  has  also  been  reported.  This  is  the  tendency 
for  an  energetic  atom  traveling  in  a  (110)  or  (100)  direction  to  be  focused 
along  a  (110)  or  (100)  line  that  is  centered  between  surrounding  (110)  or 
(100)  lines  of  lattice  atoms.  The  energetic  atom  may  travel  a  large  dis¬ 
tance  along  a  channel  compared  to  the  distance  it  would  travel  in  some 
other  direction.  The  rate  of  energy  loss  along  the  channel  is  highly 
dependent  on  the  potential  function  which  governs  the  interaction ,  the 
energy  of  the  atom,  and  its  deviation  from  the  line  of  focusing.  From 
geometrical  considerations,  one  would  expect  (100)  channels  to  dissipate 
more  energy  per  unit  length  than  would  be  dissipated  by  (110)  channels. 
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Channels,  energy  chains,  and  the  displacement  chains,  then,  may 
produce  vacancies  which  are  separated  by  large  distances  from  their 
corresponding  interstitials  and  sputtered  atoms.  Thus,  chains  and 
channels  are  competing  mechanisms  which  produce  the  same  effects . 
Energy  and  displacement  chains  should  not  be  significantly  affected  by 
the  presence  of  vacancies.  Vacancies  in  the  lines  of  atoms  surrounding 
a  channel  should  tend  to  divert  channeled  energetic  atoms  into  the 
corresponding  chains.  Interstitial  atoms  should  tend  to  disrupt  channels 
and  chains . 

It  appears  that  the  chain  and  channel  mechanisms  are  the  dominant 
factors  in  the  radiation  damage  event.  The  conditions  for  stability  and 
the  lengths  of  the  chain/channels  will  then  largely  determine  the  final 
lattice  configuration. 
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3.  Theoretical  Work  by  Others 

17 

Gibson,  et  al„  have  studied  radiation  damage  events  by  inte¬ 
gration  of  the  classical  equations  of  motion  of  the  atoms  in  a  small 
crystallite  (<C1000  atoms) .  A  central  repulsive  potential  of  the  Born- 
Mayer  type,  A  exp  (-r/b),  was  assumed  for  the  interactions  between 
atoms.  In  order  to  contain  the  events  it  was  necessary  to  restrict  the 
energy  to  <400  ev.  Calculations  were  made  on  a  high  speed  digital 
computer. 

17 

The  computer  program  developed  by  Gibson,  et  al.  keeps  a 
record  of  position  co-ordinates  and  velocity  components  for  each  atom 
in  the  crystallite.  The  orbits  of  the  atoms  are  arbitrarily  divided  into 
segments  with  respect  to  time,  i.e,  time  steps.  In  order  to  calculate 
the  trajectories  of  the  atoms  for  a  particular  time  step,  the  atoms  are 
considered  sequentially.  The  resultant  force  on  an  atom  under  con¬ 
sideration  is  calculated  from  the  position  of  the  atom  and  the  positions 
of  the  surrounding  atoms .  The  atom  is  then  moved  for  one  time  step 
under  the  influence  of  this  resultant  force.  When  the  calculations  are 
completed  for  all  atoms  for  a  given  time  step,  calculations  for  the  next 
time  step  are  begun. 

If  the  repulsive  force  between  adjacent  atoms  were  the  only  force 
simulated,  the  crystallite  would  expand  without  limit.  Therefore,  it  is 
necessary  that  an  attractive  force  act  on  the  surface  atoms  of  the  crystal¬ 
lite.  If  all  collisions  were  elastic  collisions,  and  if  the  original  energy 
(for  a  100  ev.  run)  were  completely  contained  in  the  lattice,  the  mean 
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energy  of  the  atoms  would  be  about  0.1  ev.  after  thermalization . 

This  corresponds  to  about  900°  C.  as  the  final  temperature  of  a  sys¬ 
tem  which  started  at  0°  C.  Therefore,  a  viscous  damping  force  was 
applied  to  surface  atoms  to  dissipate  energy. 

This  model  is  useful  in  studying  low  and  moderate  energy  (<400  ev.) 
effects  in  radiation  damage  events .  The  effects  studied  were  stability 
of  lattice  defects ,  anisotropicity  of  energy  transfer  (energy  chains) , 
and  displacement  chains.  Interatomic  potential  and  surface  effects 
were  also  studied. 

In  the  Gibson  model,  energy  chains  play  a  dominant  role  in  the 
dissipation  of  energy  from  the  initial  knock-on.  Frenkel  pairs  are 
shown  to  be  stable  for  distances  of  separation  greater  than  2  1/2-4 
times  the  nearest  neighbor  distance,  depending  upon  the  orientation 
of  the  pair  with  respect  to  the  lattice.  Interstitials  were  studied  ex¬ 
tensively  and  found  to  be  stable  only  in  the  split  interstitial  configura¬ 
tion  . 

21 

Oen,  et  al.  have  investigated  the  stopping  of  energetic  atoms 
(1-100  kev.)  by  solids.  The  model  is  based  on  binary  elastic  collisions 
between  a  primary  knock-on  and  successive  target  atoms  which  are 
initially  at  rest.  Lattice  effects  are  neglected;  the  target  atoms  are 
initially  randomly  distributed  throughout  the  target  material. 

The  ultimate  goal  of  the  work  by  Oen,  et  al.  is  to  discover  the 
appropriate  interatomic  potential.  The  immediate  aim  is  to  calculate 
histories  of  many  primaries  which  may  be  considered  to  simulate  a 
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beam  of  incident  atoms  or  ions. 


Trajectories  and  histories  are  calculated  only  for  the  primary 
knock-ons ,  i.e.  any  effect  produced  by  sedondary  knock-ons  is 
neglected.  The  history  of  each  primary  consists  of: 

1)  R,  the  distance  from  the  initial  position  to  the  final 
position  of  the  primary, 

2)  X,  the  component  of  R  parallel  to  the  initial  direction  of 
motion  of  the  primary, 

3)  P,  the  component  of  R  in  a  plane  perpendicular  to  the 
initial  direction  of  motion  of  the  primary,  and 

4)  L,  the  length  of  the  path  traversed  by  the  primary. 

Each  run  consists  of  the  calculation  of  histories  for  enough  primaries, 

usually  1000,  to  insure  that  statistical  fluctuations  in  the  mean  values 

of  R,  X,  P,  and  L  were  acceptable  (<  3%). 

The  potential  used  was  an  exponentially  screened  Coulomb,  or 

Bohr  potential:  V<r)  =  e*P 

*  &b 

where  and  ab=  K  4h/(2,$  -1“  7^^ 

where  a^  =  first  Bohr  radius  (0.529  A) .  The  parameter  k  may  be  adjusted 
to  achieve  agreement  with  experimental  results.  The  hard  sphere  approxi¬ 
mation  to  the  Bohr  potential  was  also  used  in  some  cases  for  comparison. 
Results  for  the  cases  of  10.5  and  60  kev.  Na^  ions  incident  upon 

aluminum  targets  were  compared  with  the  experimental  results  of  Davies  & 

0  _ 

Sims.  No  serious  discrepancy  is  found  in  the  mean  penetration ,  X . 
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However,  the  experimental  distribution  of  penetration  is  more  "skewed” 
toward  lower  penetration  than  the  distribution  of  penetration  calculated 
by  Oen,  et  al. 
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4.  The  Problem 


The  binary  collision  model  has  been  the  basis  of  all  analytical 
studies  of  radiation  damage  and  sputtering.  Many  machine  calculations 
have  been  made  in  recent  years  using  an  N  -  body  collision  model. 
However,  no  serious  attempt  has  been  made  to  determine  whether  or 
not  a  two-body  approximation  would  yield  satisfactory  results . 

Our  problem,  then,  was  to  build  a  model  based  on  the  binary  col¬ 
lision  whose  logic  would  be  sophisticated  enough  to  adequately  repre¬ 
sent  the  radiation  damage  problem.  The  nature  of  the  results  obtained 
would  then  definitely  establish  the  degree  of  confidence  to  be  placed 
in  any  binary  collision  approximation  models. 
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5.  The  Binary  Collision 


With  minor  modifications,  any  repulsive  potential  may  be  used  in 
the  two  body  interaction.  An  eroded  form  of  the  potential  is  used,  i.e. 
the  potential  has  a  definite  radius  of  effect.  We  have  taken  this  radius 
of  effect  (or  sphere  of  influence)  as  an  input  parameter,  SPHI .  There  is, 
then,  no  interaction  between  atoms  whose  separation  is  greater  than 
SPHI. 

All  applications  of  the  program  to  date  have  used  an  eroded  Born- 
Mayer  potential,  V  =  A  exp  (B (l  — r /SPHI))  for  r<  SPHI.  The  constants 
A  &  B  were  taken  from  Gibson  et  al.  for  three  different  Born-Mayer 
potentials,  termed  simply  Potentials,  1,2,  and  3.  The  constants  for 
these  three  potentials  are: 

POTENTIAL  A  (ev.)  B 

1  .0392  16.97 

2  .0510  13.00 

3  .1004  10.34 

The  interaction  starts  with  the  two  atoms  at  SPHI  distance  of  sepa¬ 
ration  with  all  the  energy  as  kinetic  energy.  All  co-ordinates  and  velo¬ 
city  components  are  transformed  to  the  center  of  mass  system,  where 
the  equations  of  motion  are: 

1)  dt=(2/m'  (E-V-L2/2m'r2))-1/2  dr 

2)  d0  =  (L/m'r2)  dt 
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where  t  =  time 

m'  =  reduced  mass 
E  =  kinetic  energy 
V  =  potential  energy 
L  =  total  angular  momentum 
r  =  distance  of  separation 
0  =  angle  of  deviation 
s  =  impact  parameter 

The  equations  are  integrated  from  r  =  SPH1  to  r  =  CPA,  the  distance 
of  closest  approach.  Since  the  interaction  is  symmetrical  in  the  center 
of  mass  system,  the  definite  integrals  of  the  above  equations  are  one- 
half  the  time  of  interaction  and  one-half  the  total  angular  deviation, 
respectively . 

22 

CPA  is  found  by  an  interation  of  the  formula: 

(CPA)2  (1  -  V(CPA/Er)  -  S2  =  0 

The  first  approximation  to  CPA  is  taken  as  SPHI.  A  four  point  Gauss- 

23 

Legendre  quadrature  is  employed  to  evaluate  the  integrals. 
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6  .  Development  of  the  Model 


When  the  possibility  of  a  feasible  two  body  approximation  to  the 
radiation  damage  problem  is  considered,  two  major  problems  are  im¬ 
mediately  encountered . 

The  atoms  would,  of  necessity,  be  considered  in  some  sort  of 
ordered  sequence,  but  all  calculations  for  a  single  atom  could  not  be 
carried  through  before  the  other  atoms  are  considered  since  this  would 
obviously  leave  out  of  the  calculations  the  effect  of  collisions  suffered 
by  other  atoms.  This  first  major  problem  led  to  the  time  step  approach, 
and  eventually  to  the  use  of  an  absolute  time  of  an  atom's  last  collision. 

The  second  major  problem  is  to  develop  a  mechanism  which  will 
handle  simultaneous  collisions  within  a  two-body  approximation. 

After  the  program  was  developed,  we  discovered  that  the  lattice 
would  not  hold  together  which  necessitated  a  further  modification. 

This  difficulty  is  only  partially  resolved. 

The  first  question  considered  in  the  development  of  the  model  was 
to  mathematically  create  a  face  centered  cubic  lattice  in  table  form, 
together  with  a  method  which  tags  each  atom  so  that  its  subsequent 
movements  and  initial  position  could  be  found.  Figure  2,  APPENDIX  II, 
shows  the  ordered  arrangement  found  in  a  (100)  plane  of  a  fee  lattice. 

The  distance  between  two  atoms  along  the  (100)  direction  is  taken  as 
two  units.  It  will  be  noted  that  if  the  total  X  dimension,  IX,  is  odd, 
then  there  are  an  equal  number  of  atoms  in  the  line  y  =  0  and  the  line 
y  =  1 .  By  specifying  EX  odd ,  the  number  of  atoms  in  a  line  y  =  constant 
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is  known  to  be  (EX  +  l)/2.  This  value,  NLINE,  can  then  be  multiplied 
by  the  total  Y  dimension  plus  one  (IY+1) ,  and  the  number  of  atoms  in 
the  whole  plane,  NPLANE,  is  obtained.  It  will  be  noted  from  Figure  3 
that  the  next  plane,  one  Z  unit  removed,  is  akin  to  a  photographic 
negative  of  the  plane  under  consideration.  However,  the  ideas  of 
NLINE  and  NPLANE  are  still  valid .  The  total  number  of  atoms  in  the 
lattice  can  then  be  found  by  multiplying  NPLANE  by  the  Z  dimension 
plus  one  (IZ+1).  If  the  atoms  are  numbered  sequentially  as  shown  in 
APPENDIX  II,  the  initial  co-ordinates  of  their  sites  can  also  be  found 
by  the  method  shown  in  APPENDIX  II  and  APPENDICES  IV  &  V  (BOX  31) . 

The  process  outlined  above  will  create  and  number  a  1500  atom 
lattice  in  approximately  one  second.  An  advantage  of  this  method  is 
that  a  rectangular  lattice  of  any  size  or  shape  desired  can  be  created 
by  specifying  IX,  IY,  and  IZ.  The  only  restriction  is  that  IX  must  be 
odd . 

The  atoms  are  now  considered  sequentially,  for  a  period  of  one 
Time  Step  Length  (TSL) ,  and  all  collisions  that  will  take  place  before 
the  end  of  the  time  step  are  allowed  to  happen .  All  atoms  are  then 
moved  to  points  in  space  and  time  corresponding  to  the  end  of  the  time 
step,  and  the  process  is  repeated  for  the  next  time  step. 

The  length  of  a  time  step  was  originally  calculated  so  that  the 
most  energetic  atom  in  the  lattice  could  move  one  lattice  unit  (1.807 
Angstroms  for  copper)  in  one  time  step.  This  has  been  modified  so 
that  it  is  now  able  to  move  any  desired  fraction,  TFAC ,  of  one  lattice 
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unit.  TFAC  is,  naturally,  one  of  the  input  parameters  of  the  problem. 

Since  the  energy  of  the  most  energetic  particle,  ENERGY,  in  the 
lattice  decreases  during  the  course  of  the  run,  the  TSL  should  be  re¬ 
calculated  periodically  for  greater  efficiency.  Provisions  were  made 
to  recalculate  the  Time  Step  Length  every  JPB  (an  input  parameter)  time 
steps.  The  TSL  is  also  recalculated  every  time  another  particle  is  shot 
into  the  lattice,  since  it  is  assumed  that  this  particle  will  be  the  most 
energetic  one  present,  and  it  was  not  thought  feasible  to  wait  for  the 
next  regular  recalculation.  A  table  of  the  first  30  Time  Step  Lengths, 
and  the  time  step  of  recalculation,  is  kept  for  reference. 

As  now  constituted,  there  are  one  external  and  two  internal  methods 
used  to  determine  the  number  of  time  steps  which  complete  a  run.  The 
external  method  is  the  only  one  available  to  the  operator  once  a  com¬ 
putational  run  has  commenced. 

First,  a  variable  limit  is  set  on  the  maximum  number  of  time  steps 
(MNTS  ,  an  input  parameter)  that  may  be  taken .  The  calculation  must 
cease  when  this  number  of  time  steps  has  been  completed,  and  the 
results  are  then  obtained . 

Second,  the  program,  at  the  beginning  of  each  time  step,  decides 
whether  or  not  it  should  continue  for  another  time  step.  The  decision 
to  continue  is  made  if  the  energy  of  the  most  energetic  particle  (called 
ENERGY)  in  the  lattice  (excluding  all  those  that  have  left  the  lattice) 
is  greater  than  a  certain  threshold  energy,  THERMAL  (an  input  parameter 
explained  below) .  If  ENERGY  is  lesli  than  THERMAL,  the  run  is 
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considered  completed  and  the  informational  output  is  given. 

Third ,  the  operators  may  decide  to  shut  down  operation  before 
either  of  the  two  internal  limits  have  been  reached.  This  may  be  done 
at  any  time,  simply  by  throwing  a  switch,  and  the  computation  will 
cease  at  the  beginning  of  the  next  time  step. 

The  external  mode  of  problem  completion  should  not  be  the  normal 
operating  condition,  nor  should  the  number  of  time  steps  be  allowed  to 
approach  MNTS  under  most  conditions,  since  either  method  still  leaves 
atoms  in  the  lattice  with  more  energy  than  the  quantity  THERMAL. 
Usually  the  calculation  should  cease  of  its  own  accord,  with  the 
second  method  given  above  as  the  basis  for  the  internal  decision. 

Other  factors,  such  as  running  time,  mechanical  failures,  etc.  ,  may 
necessitate  the  use  of  the  external  method,  but  these  are  to  be  avoided 
if  possible . 

A  method  of  program  regeneration  was  also  devised  so  that  the  run 
could  be  discontinued  at  the  start  of  any  time  step  and  then  restarted  at 
that  point  at  some  later  time.  This  has  definite  value  where  continuous 
periods  of  available  machine  time  are  not  long,  or  when  mechanical  or 
electrical  failures  destroy  the  calculations .  The  information  needed  for 
regeneration  can  be  obtained  at  the  beginning  of  any  time  step  and  the 
calculation  continued  (or  stopped  if  desired) . 

In  this  manner,  if  the  calculation  has  been  running  for  three  hours 
before  a  failure  of  some  kind  occurs ,  the  only  time  that  need  be  re¬ 
peated  on  the  machine  is  the  time  elapsed  since  the  last  regeneration 
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information  was  put  on  magnetic  tape.  It  should  be  noted  that  the  re¬ 


generation  need  not  be  undertaken  at  once ,  but  may  be  delayed  until 
convenient.  This  regeneration  process  has  great  practical  importance 
to  the  user  of  the  program,  but  is  of  no  theoretical  value. 

Many  atoms  are  given  only  a  cursory  examination  during  a  time 
step,  while  the  more  energetic  atoms  are  subjected  to  a  highly  detailed 
treatment.  The  use  of  two  energy  thresholds  in  the  model  simplifies 
this  choice  for  a  particular  atom. 

ETHRESH,  the  upper  energy  threshold,  is  the  energy  required  for 

an  atom  to  move  from  its  site  to  a  site  near  it,  i.e.  the  displacement 

energy.  Calculations  by  others ^  have  shown  this  displacement 

threshold  to  exist,  and  also  that  it  is  in  the  neighborhood  of  20 

1 7 

electron  volts.  They  have  also  shown  that  this  displacement  thresh¬ 
old  depends  upon  direction,  and  is  lower  in  the  (100)  and  (110)  direc¬ 
tions  than  in  other  directions  in  the  lattice.  This  directional  depend¬ 
ence  of  the  displacement  threshold  has  not  been  built  into  the  present 
model,  but  it  should  be  considered  in  any  future  development  of  the 
model . 

THERMAL,  the  lower  energy  threshold,  is  an  arbitrary  dividing 
line  imposed  upon  the  model,  below  which  atoms  are  quite  arbitrarily 
placed  on  sites  or  formed  into  interstitial  pairs  with  other  atoms  and 
not  allowed  to  move.  Atoms  above  THERMAL  are  allowed  to  participate 
normally  in  collisions . 
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There  are,  then,  three  distinct  energy  classes  among  the  atoms  in 


the  lattice: 

1)  Those  with  energy  above  ETHRESH, 

2)  Those  with  energy  between  ETHRESH  and  THERMAL,  and 

3)  Those  with  energy  below  THERMAL. 

The  three  classes  will  be  discussed  separately. 

The  atoms  above  ETHRESH  do  not  occupy  a  site  in  the  lattice, 
collide  with  other  atoms  as  necessary,  and  obtain  information  about 
their  subsequent  motion  only  from  the  interaction  subroutine . 

The  atoms  between  ETHRESH  and  THERMAL  are  forced  to  occupy 
specific  lattice  sites,  but  are  allowed  to  vibrate  around  these  sites. 
They  cannot  move  to  a  new  site  because  their  energy  is  less  than  the 
displacement  energy.  These  atoms  are  allowed  to  move  undisturbed 
until  they  have  a  collision,  and  then  after  the  collision  has  been  com¬ 
pleted,  if  they  lose  a  specific  percentage,  TURN  (an  input  parameter), 
of  their  original  energy  their  velocity  vectors  are  changed  to  point  back 
toward  the  site  occupied  by  the  atom. 

We  found  prior  to  the  use  of  this  turn-around  method  that  the  lattice 
had  no  inherent  cohesiveness.  This  method  holds  the  lattice  together. 

To  speed  up  computation  time,  atoms  with  energy  below  THERMAL 
are  forced  to  occupy  a  fixed  site,  or  to  share  an  interstitial  site  with  a 
partner,  and  are  not  allowed  to  move.  This  is  a  reasonable  and  valuable 
approximation  for  atoms  of  low  energy,  since  they  will  not  have  an 
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appreciable  effect  on  the  results  of  the  computation.  The  energy  and 
associated  velocity  components  of  these  thermalized  atoms  are  kept 
intact  for  use  in  possible  future  collisions. 

After  the  decision  is  made  to  continue  for  another  time  step,  and  the 
Time  Step  Length  has  been  recalculated,  the  atoms  are  considered  sequen¬ 
tially,  The  information  associated  with  the  atom  under  consideration 
(primary  atom)  which  is  stored  in  the  LB  array  is  first  inspected.  If  the 
atom  has: 

*  •  •  i 

l 

1)  An  energy  equal  to  or  less  than  THERMAL,  or  has 

2)  Been  through  the  current  time  step ,  or  has 

3)  Left  the  lattice , 

we  move  on  to  consideration  of  the  next  atom. 

We  will  account  for  all  collisions  if  we  consider  only  those  atoms 
with  energy  greater  than  THERMAL,  since  atoms  with  energy  less  than 
THERMAL  are  not  allowed  to  move  unless  hit.  Atoms  that  have  already 
been  through  the  time  step  have  completed  all  interactions  which  occur 
during  the  time  step  and  should  not  become  involved  again  until  the 
next  time  step.  Atoms  that  have  left  the  lattice  will  not  be  considered 
for  the  remainder  of  the  problem . 

A  primary  atom  with  energy  above  THERMAL  is  finally  selected,  or 
the  computation  shuts  down.  A  list  of  all  other  atoms  it  is  possible  for 
this  atom  to  hit  in  this  time  step  is  required.  To  shorten  the  computa¬ 
tion  ,  a  mathematical  box  is  constructed  around  the  atom  in  question 
(hereafter  referred  to  as  atom  M)  and  a  list  (the  NUM  list)  is  made  of 
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all  atoms  whose  positions  (centers)  are  inside  the  box.  If  the  box  is 
made  large  enough,  it  will  be  impossible  for  atom  M  to  hit  any  atoms 
outside  the  box,  so  we  need  consider  only  those  on  the  NUM  list.  At 
present,  the  box  is  given  sides  of  3.02  units,  orjf  1.51  units.  Initially, 
the  sides  were  +  1.01,  but  this  has  been  found  to  be  too  small,  and 
obvious  collisions  were  missed.  Atoms  that  have  been  through  the  time 
step  are  not  included  on  this  list,  since  their  co-ordinates  are  relative 
to  the  end  of  the  time  step  rather  than  the  start  of  the  time  step.  Atom  M 
is  placed  first  on  the  list  for  programming  facility  and  easy  reference. 

The  impact  parameter  between  every  other  atom  on  the  NUM  list  and 
atom  M  is  then  calculated,  these  distances  comprise  the  DSTANCE  list. 
The  time  at  which  each  atom  reaches  this  distance  of  closest  approach, 
and  the  time  at  which  the  atom  and  atom  M  are  SPHI  (sphere  of  influence 
of  an  atom)  distance  apart  are  also  calculated.  These  times  are  placed 
on  the  TMIN  and  T  lists,  respectively. 

During  the  course  of  the  calculation  of  the  DSTANCE,  TMIN,  and  T 
lists,  several  tests  are  performed  on  these  quantities.  If  the  DSTANCE 
for  an  atom  is  larger  than  SPHI,  no  collision  will  occur.  If  the  relative 
velocity  is  less  than  10“'®,  then  the  atoms  are  moving  very  slowly 
relative  to  one  another,  and  the  collision  is  ignored.  Several  tests  are 
then  performed  on  each  T.  If  the  time  T  is  less  than  the  start  of  the  pre- 
ceeding  time  step,  it  would  have  been  considered  in  this  previous  time 
step,  and  is  therefore  treated  as  a  spurious  collision.  If  the  time  T  is 
greater  than  the  end  of  the  current  time  step,  it  will  be  re-examined  in 
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the  next  time  step.  The  times  of  the  last  collision  of  atom  M  and  the 
other  atom  are  then  checked,  and  if  the  time  T,  on  an  absolute  scale, 
is  less  than  either  of  these,  the  collision  is  spurious  since  it  would  be 
starting  before  one  that  had  already  taken  place.  The  last  test  performed 
is  whether  the  time  T  is  within  .00001  (TSL)  of  the  end  of  the  time  step. 
If  so,  we  feel  that  it  would  be  better  to  consider  this  collision  in  the 
next  time  step,  since  to  consider  it  now  would  mean  that  the  end  of  the 
collision  was  far  into  the  next  time  step,  which  is  unsatisfactory  from 
a  practical  viewpoint. 

If,  after  this  series  of  eliminations,  there  are  no  atoms  left,  then 
atom  M  is  moved  to  the  end  of  the  time  step,  and  the  next  atom  in  se¬ 
quence  becomes  atom  M. 

If  atoms  remain  on  the  list  we  select  the  ones  that  atom  M  will  hit. 
Two  alternative  methods  of  selection  were  considered  by  the  authors: 

1)  Select  the  one,  or  ones,  that  are  closest  (i.e.  a  space- 
wise  selection) ,  or 

2)  Select  those  that  will  collide  with  atom  M  first  (i.e.  a 
time-wise  selection) . 

Only  the  second  method  has  been  used. 

The  space-wise  selection  would  surely  include  all  the  hardest 
collisions,  and  would  probably  be  easier  from  a  programming  viewpoint. 
However,  a  time-wise  selection  insures  that  atom  M  will  hit  first  what 
it  really  should  hit  first.  Atom  M  may  thereby  miss  a  hard  collision 
with  another  atom  because  of  another  previous  soft  one ,  but  we 
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currently  feel  that  this  is  a  truer  picture  of  what  actually  happens  in  the 
lattice  during  the  radiation  damage  process.  We  feel  that  both  methods 
should  be  tried  in  future  development  of  the  model.  A  study  as  to  the 
relative  feasibility  of  both  methods  would  enhance  the  usefulness  of  the 
model . 

Atoms  are  selected  that  have  the  smallest  times  ,  T.  Since  we  assume 
that  two  atoms  have  no  interaction  if  they  are  separated  by  a  distance 
greater  than  SPHI,  the  time  T  is  the  actual  time  of  collision.  We  feel 
that  some  leeway  should  be  available  here ,  therefore  all  atoms  whose 
times  T  are  within  a  certain  amount  (CUTOFF)  of  a  time  step  from  the 
minimum  times  are  also  selected.  These  atoms  are  placed  on  the  KHIT 
list,  with  atom  M  first  on  the  list.  Since  CUTOFF  is  an  input  parameter, 
adjustment  to  a  satisfactory  value  may  be  accomplished  easily.  Use  of 
the  CUTOFF  parameter  also  eliminates  some  of  the  objections  to  the 
time-wise  selection  process. 

The  initial  velocity  components  and  energies  of  all  atoms  on  the 
KHIT  list  are  then  stored  in  the  SAVE  array  because  they  will  be  needed 
for  the  velocity  scaling  processes. 

It  must  be  emphasized  that  atom  M  hits  all  those  on  the  KHIT  list 
almost  simultaneously  „  This  necessitates  some  sort  of  energy  and 
velocity  scaling  process,  or  a  general  solution  to  the  N  body  problem. 

We  feel  that  if  the  two  body  approximation  is  to  have  any  validity 
atom  M  must  enter  into  each  of  the  two  body  collisions  with  its  initial 
energy,  but  the  program  must  scale  the  results  in  some  fashion.  Assume, 
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for  example,  that  M  initially  has  an  energy  of  100  ev.  and  hits  two 
atoms  simultaneously  with  identical  impact  parameters.  It  is  more 
feasible  to  assume  that  M  hits  each  one  with  an  energy  of  100  ev. 
than  to  assume  that  M  hits  one  with  100  ev.  and  hits  the  other  with 
whatever  energy  is  left,  say  37  ev. 

Naturally,  with  M  hitting  each  atom  separately  with  100  ev.  the 
struck  atoms  will  receive  too  much  energy,  we  must  therefore  scale  the 
results . 

In  the  actual  collision  process  of  the  model,  then,  atom  M  and  the 
next  atom  on  the  KHIT  list  are  moved  to  their  calculated  positions  at  the 
time  T  for  the  atom  being  hit.  The  collision  is  then  calculated  and  the 
struck  atom  is  given  its  new  position  co-ordinates,  velocity  components, 
and  energy.  Atom  M  and  the  next  atom  on  the  KHIT  list  are  then  moved 
to  the  positions  corresponding  to  their  actual  time  of  collision,  T,  and  so 
forth.  This  is  repeated  until  atom  M  has  "struck"  all  the  atoms  on  the 
KHIT  list. 

To  conserve  both  energy  and  momentum  in  any  scaling  process 
would  constitute  a  general  solution  to  the  N  body  problem.  We,  there¬ 
fore,  elected  to  conserve  energy  rather  than  momentum,  feeling  that 
this  would  give  us  a  more  realistic  insight  into  the  radiation  damage 
problem. 

There  are  four  general  scaling  cases  that  must  now  be  considered: 

1)  Atom  M  "lost"  more  energy  than  it  started  with, 

2)  Atom  M  "lost"  exactly  all  the  energy  it  started  with, 
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3)  Atom  M  started  with  an  energy  above  ETHRESH  and  did  not 

"lose"  all  of  its  energy,  and 

4)  Atom  M  started  with  an  energy  below  ETHRESH,  but  above 
THERMAL,  and  did  not  lose  ail  of  its  energy. 

These  cases  will  be  discussed  in  order. 

CASE  I 

Atom  M  is  allowed  to  lose  all  its  initial  energy,  and  its  velocity 
components  and  energy  are  set  equal  to  zero.  The  energies  and 
velocities  of  the  struck  atoms  are  then  scaled  to  conserve  energy.  The 
energy  scaling  factor  used  in  this  case  is  simply  the  total  energy  (of 
those  on  the  KHIT  list)  before  the  collisions  divided  by  the  total  energy 
after  the  collisions.  The  velocity  scaling  factor  is  the  square  root  of 
the  energy  scaling  factor.  The  energies  and  velocities  of  all  atoms  on 
the  KHIT  list ,  except  atom  M ,  are  scaled  by  these  two  scaling  factors . 

Both  energy  and  the  absolute  velocity  directions  are  conserved  by 
this  process. 

CASE  II 

Atom  M  is  allowed  to  lose  all  of  its  initial  energy,  and  its  velocity 
components  and  energy  are  set  equal  to  zero  Here,  however,  energy 
scaling  is  not  necessary  since  energy  has  already  been  conserved. 

CASE  III 

Here  it  is  not  necessary  to  scale  the  energies  of  the  struck  atoms, 
since  atom  M  did  not  lose  all  of  its  initial  energy,  so  the  struck  atoms 
retain  the  energies  and  velocities  they  received  in  the  two  body  collision 
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process.  The  energy  of  M  is  found  by  subtracting  the  energy  gained  by 
the  struck  atoms  from  M's  initial  energy.  M’s  velocity  components 
must  now  be  found  , 

The  sums  of  the  individual  velocity  components  calculated  for  M, 
but  not  given  to  M,  are  taken.  Call  these  sums  S^,  S2,  and  S3  for  the  X, 
Y,  and  Z  directions.  Let  R  be  (Sj  +  The  magnitude  of  M's 

final  velocity,  V,  is  known  since  the  energy  is  known.  The  X  velocity 

I 

is  then  taken  as  (V/R)  (Sj),  and  similarly  the  Y  and  Z  velocities  are 
(V/R)^)  and  (V/RHS3)  respectively. 

CASE  IV 

Two  sub-cases  must  be  considered.  First,  if  atom  M  does  not  lose 
TURN  %  of  its  initial  energy,  the  treatment  used  in  CASE  III  will  be 
followed . 

Second,  if  atom  M  does  lose  a  significant  portion  of  its  energy,  its 
velocity  components  are  changed  so  that  the  atom  is  directed  back  toward 
the  site  it  occupies .  The  energy  of  M  is  again  its  initial  energy  minus 
the  energy  gained  by  the  struck  atoms.  The  scale  factor  used  is 
(Velocity/distance  from  site).  This  is  multiplied  by  AX,  AY,  and 
A  Z  to  give  the  velocity  components .  ^X  is  simply  the  X  co-ordinate 
of  the  site  minus  the  X  co-ordinate  of  M. 

In  all  cases,  atom  M's  spatial  co-ordinates  are  found  by  averaging 
the  individual  co-ordinates  calculated  for  M  in  the  X,  Y,  and  Z  directions 
(i.e.  dividing  the  sum  by  NHIT) . 
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If  atom  M  is  a  primary  atom  (i.e.  one  of  those  considered  sequen¬ 
tially  for  this  time  step) ,  then  the  KHIT  list  is  duplicated  by  the  LAST 
list.  Atom  M  will  then  proceed  to  go  through  its  designated  collisions, 
look  for  more  collisions,  and  go  through  these  until  it  completes  all 
collisions  in  this  time  step. 

When  atom  M  has  completed  all  its  possible  collisions  for  this  time 
step,  each  atom  that  M  hit  (found  in  order  on  the  LAST  list)  undergoes 
the  same  process  and  attempts  to  hit  other  atoms.  The  atoms  hit  by  M, 
or  others  on  the  LAST  list,  are  not  inspected  for  possible  further  colli¬ 
sions  in  this  time  step.  We  feel  that  these  "secondary"  atoms  would 
not  have  enough  time  to  have  a  collision  in  this  time  step,  and  if  they 
did,  the  collisions  will  be  calculated  in  the  next  time  step.  At  this  point, 
a  possibility  does  exist  that  some  reasonable  collisions  are  excluded 
by  the  process,  but  extensive  calculations  have  not  revealed  any 
examples  of  this  exclusion  to  date. 

The  absolute  time  of  an  atom's  last  collision  is  also  calculated 
during  the  collision  processes  described  above.  For  each  of  the  struck 
atoms  this  time  is  the  total  time  from  the  start  of  the  computation  to  the 
start  of  this  time  step  plus  the  time,  T,  of  collision  (T  is  measured  from 
the  start  of  the  current  time  step) . 

The  time  of  M's  last  collision  is  taken  as  the  start  of  its  collision 
with  the  last  atom  on  the  KHIT  list  (computed  as  above)  plus  a  certain 
percentage  (PERCENT ,  an  input  parameter)  of  the  present  time  step 
length.  M  is  thus  prevented  from  participating  in  further  collisions 
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for  a  definite  period  of  time  after  the  start  of  its  last  collision.  We 
feel  that  this  approach  is  a  valid  one ,  however  study  of  this  aspect  of 
the  model  should  be  continued  to  determine  its  true  usefulness . 

Sites  must  now  be  found  for  those  atoms  on  the  KHIT  list  whose 
pre-collision  energies  were  above  ETHRESH  and  whose  energies  are 
now  below  ETHRESH. 

If  the  rounded-off  co-ordinates  of  an  atom  in  this  class  (call  it  J) 
give  the  co-ordinates  of  an  actual  lattice  site ,  then  the  atom  will 
occupy  that  site,  possibly  as  part  of  an  interstitial  pair  if  the  site  is 
already  occupied.  Atom  J  will  then  be  allowed  to  vibrate  around  its 
site  until  its  energy  falls  below  the  THERMAL  cutoff.  If  J's  rounded-off 
co-ordinates  are  not  the  co-ordinates  of  a  site,  then  the  nearest  neighbor 
sites  are  found  and  the  site  closest,  space-wise,  to  atom  J  is  chosen  as 
J's  vibrational  center.  Again,  this  may  result  in  the  formation  of  an 
interstitial  pair.  Since  one  of  the  atoms  has  an  appreciable  energy, 
however,  this  interstitial  pair  should  be  destroyed  in  the  next  time  step. 

The  site  decided  upon  may  already  be  occupied  by  an  interstitial 
pair:  if  so,  a  "triple  interstitial"  is  formed.  An  error  counter  is  then 
increased  (MISTAKE  (5)).  The  possibility  of  "triple  interstitials"  could 
be  eliminated  by  more  extensive  programming,  but  an  evaluation  should 
first  be  made  of  the  seriousness  of  this  error.  The  frequency  of  forma¬ 
tion  of  these  "triple  interstitials"  has  been  small  thus  far,  not  more 
than  two  or  three  per  run,  which  is  tolerable.  If  this  frequency  reaches 
an  intolerable  level  in  future  development,  the  model  must  be  modified 
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to  eliminate  any  possibility  of  their  occurrence. 


All  atoms  on  the  KHIT  list  whose  post-collision  and  scaling  energies 
are  below  THERMAL  are  placed  on  the  LATER  list  at  this  time.  Use  of 
the  LATER  list  is  discussed  below. 

If  one  of  the  struck  atoms  on  the  KHIT  list  was  a  member  of  an 
interstitial  pair,  the  other  atom  must  now  be  considered.  If  this  other 
atom  was  not  hit,  and  the  struck  atom  now  has  an  energy  above  ETHRESH 
(a  very  unlikely  occurrence) ,  then  the  other  atom  is  placed  back  on  the 
site,  rather  than  being  left  in  the  split  interstitial  position.  We  feel 
that  moving  the  atom  back  onto  the  site  is  more  realistic  than  the 
definite  error  of  leaving  it  in  the  split  interstitial  position. 

All  calculations  in  the  model  are  made  on  the  assumption  that  the 
atoms  are  both  space-wise  and  time-wise  at  the  start  of  the  current  time 
step,  except  those  that  have  completed  the  current  time  step.  Therefore, 
all  the  atoms  on  the  KHIT  list  are  now  moved  back  in  time  and  space  on 
their  post-collision  tracks  to  the  start  of  this  time  step . 

The  entire  collision  process,  including  the  construction  of  a  new 
NUM  list  is  now  repeated  for  atom  M.  This  repetition  will  continue  for 
M  until  M  is  unable,  time-wise,  to  have  further  collisions  in  this  time 
step.  The  next  atom  on  the  LAST  list  will  then  become  atom  M  and  the 
collision  process  will  be  repeated  until  this  M  is  unable  to  have  further 
collisions  in  this  time  step.  This  overall  procedure  is  then  repeated 
until  the  LAST  list  has  been  exhausted. 
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All  the  atoms  on  the  LAST  list,  except  those  on  the  LATER  list, 

(the  original  M  is  the  first  atom  on  the  LAST  list)  are  now  moved,  space 
and  time-wise,  to  the  end  of  the  time  step.  These  atoms  have  now  com¬ 
pleted  the  time  step  and  will  not  enter  into  the  subsequent  calculations 
until  the  next  time  step . 

The  LATER  list  is  now  utilized  to  place  the  atoms  that  have  partici¬ 
pated  in  this  sequence  of  collisions  and  whose  energies  are  below 
THERMAL  on  sites,  either  as  members  of  interstitial  pairs  or  by  them¬ 
selves  on  previously  vacant  sites. 

The  co-ordinates  of  the  atoms  on  the  LAST  list  are  now  examined . 

If  an  atom  has  left  the  lattice,  the  time  step  number  and  the  particular 
face  (see  APPENDIX  I,  NFACE)  of  exit  is  recorded  in  the  LB  array. 

This  completes  the  consideration  of  atom  M  (the  original  M)  for  this 
time  step.  The  next  atom  (M  +  1)  is  then  considered  in  the  same  manner 
until  all  atoms  have  been  through  the  procedure  (and,  consequently,  the 
time  step) . 

After  all  the  atoms  have  been  through  these  processes ,  all  un¬ 
affected  atoms  are  moved  to  the  end  of  the  time  step,  which  completes 
the  time  step.  A  new  time  step  is  then  started,  or  the  run  is  shut  down 
by  one  of  the  three  processes  mentioned  earlier. 

A  numerical  example  of  the  order  in  which  atoms  are  considered  may 
be  of  some  value  here.  Assume  atom  #73  (M)  hits  four  atoms  simulta¬ 
neously,  atoms  75,  79,  87,  and  65.  Assume  further  that  #73  hits  them 
in  the  order:  #79,  #65,  #87,  #75.  Atom  #73  will  hit  these  four,  then 
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repeat  the  process  for  further  collisions.  Assume  #73  hits  #80  and  #67 
next,  but  has  no  further  collisions  after  that.  Then,  atom  #79  will 
become  the  next  M  since  it  is  second  on  the  LAST  list,  #65  the  next  M, 
and  after  #75  has  completed  its  possible  collisions,  if  any,  the  regular 
sequence  will  be  taken  up  once  more.  Atom  #73  was  the  last  one  con¬ 
sidered  sequentially,  therefore  atom  #74  now  becomes  the  primary  atom, 
atom  M . 

As  mentioned  above ,  atoms  whose  energies  fall  below  THERMAL  are 
either  placed  alone  on  vacant  sites ,  or  formed  into  split  interstitial  pairs 
with  other  atoms  near  them.  This  is  done  only  after  a  primary  atom,  to¬ 
gether  with  those  on  its  LAST  list,  has  been  completely  considered.  If 
this  were  done  at  the  end  of  every  collision  process,  it  would  severely 
hamper  the  multiple  collision  process,  and  Impose  a  great  deal  of  strain 
on  the  credibility  of  the  argument  behind  it. 

When  a  list  of  the  nearest  neighbor  sites  is  formed,  if  the  atom  in 
question  (call  it  MJ)  is  within  one  unit  of  a  face  of  the  crystallite ,  some 
of  the  geometrically  nearest  neighbors  normally  found  will  not  be  present. 
There  are  22  special  cases  for  each  one  of  the  two  general  cases  (MJ's 
co-ordinates  are  a  site  or  non-site) .  Sixteen  of  these  special  cases 
for  each  concern  the  edges  and  comers  of  the  crystal.  In  our  model, 
the  32  (total)  possible  cases  for  atom  MJ  near  an  edge  or  comer  are  not 
considered  separately,  but  are  Included  in  the  12  special  cases  of  when 
MJ  is  near  a  lattice  face.  When  one  of  the  32  unconsidered  special 
cases  does  arise,  an  error  counter  is  increased.  Inspection  of  all 
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preliminary  results  to  date  have,  however,  shown  NO  discrepancies  as 
a  result  of  this  simplification  at  the  edges  and  corners  of  the  lattice. 

After  the  list  of  nearest  neighbor  sites  ,  the  MOON  list,  is  formed, 
these  sites  are  inspected  for  possible  vacancies.  If  only  one  vacancy 
exists,  atom  MJ  is  placed  on  that  site  and  the  appropriate  corrections 
are  made  to  the  LB  array.  If  more  than  one  vacancy  exists  among  MJ's 
nearest  neighbor  positions,  then  one  of  these  must  be  selected  for  MJ's 
occupancy.  We  feel  that  this  selection  should  be  made  on  a  space-wise 
condition  of  proximity,  rather  than  on  any  time-wise  condition.  There¬ 
fore,  atom  MJ  is  allowed  to  occupy  the  closest  vacant  site.  If  there  is 
more  than  one  vacancy  at  this  smallest  distance,  MJ  will  be  placed  on 
its  own  site,  if  available,  or  the  selection  will  be  made  in  a  random 
manner.  The  choice  of  a  particular  vacancy  in  the  case  of  multiple 
vacancies  is  somewhat  unsophisticated,  however  at  this  stage  of  devel¬ 
opment  in  the  radiation  damage  field,  we  feel  that  no  process  of  a  more 
sophisticated  nature  can  be  justified. 

If  no  vacancies  exist  among  atom  MJ's  nearest  neighbor  sites,  then 
a  prospective  interstitial  site  and  partner  are  chosen.  If  atom  MJ  has 
no  energy  or  velocity,  this  decision  is  made  on  a  space-wise  basis, 
i.e.  the  closest  site  is  chosen.  If  atom  MJ  does  have  energy,  then  the 
selection  is  made  on  a  time-wise  basis,  i.e.  the  site  that  atom  MJ 
passes  first. 

Once  the  selection  of  a  prospective  interstitial  partner  has  been 
completed,  the  atom  (call  it  MM)  which  occupies  the  site  is  found. 
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Atom  MM's  nearest  neighbor  sites  are  then  found,  and  the  list  is 
examined  for  vacancies.  If  any  vacancies  exist,  then  atom  MM  is 
moved  to  one  of  these  vacant  sites  by  the  same  method  discussed 
above  for  MJ.  Atom  MJ  is  then  moved  onto  atom  MM’s  recently  vacated 
site  and  no  interstitial  is  formed. 

If,  however,  atom  MM  finds  no  vacancies  ,  then  an  interstitial  pair 

must  be  formed.  We  feel  that  of  the  three  types  of  interstitial  pairs 

mentioned  earlier,  the  split  interstitial  is  the  most  plausible.  Indeed, 

Gibson,  et  al.  have  shown  the  split  interstitial  to  be  the  only  stable 
1 7 

configuration.  The  AX,  AY,  and  £Z  distances  between  atom  MJ 
and  the  site  atom  MM  occupies  are  found.  The  largest  of  these  is 
taken  as  the  axis  of  the  interstitial  pair.  Each  atom  is  then  placed  on 
its  split  position,  0.5  units  from  the  site  along  the  axis.  Once  the 
appropriate  corrections  are  made  to  the  LB  array,  the  formation  of  the 
pair  is  completed. 
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7 .  Operating  Procedures 


The  model  consists  of  three  separate  but  inter-connected  programs, 
MASTER,  SLAVE,  and  RON.  Machine  memory  limitations  forced  the  devel¬ 
opment  of  three ,  rather  than  one ,  programs . 

The  main  program,  MASTER,  is  the  computational  program  and  does 
all  of  the  calculations  necessary  to  complete  the  run.  The  two  second¬ 
ary  programs,  SLAVE  and  RON,  merely  use  the  binary  output  of  program 
MASTER  as  input  and  put  the  results  in  a  semi-analyzed  readily  legible 
form.  The  orbits  (or  intermediate  positions  of  the  moving  atoms)  are 
given  by  program  RON  while  program  SLAVE  tabulates  and  categorizes 
the  final  positions,  energies,  etc.  of  the  lattice  atoms.  Again,  be¬ 
cause  of  memory  limitations,  the  two  secondary  programs  may  not  be 
combined . 

If  a  binary  version  of  program  MASTER  is  obtained  with  the  FORTBIN 
compiler,  the  storage  required  may  be  found  by  listing  of  the  binary  tape. 
The  maximum  lattice  may  then  be  computed .  The  FORTMAP  compiler 
should  be  used  to  compile  the  large  version  of  MASTER  since  it  requires 
less  memory  than  the  other  compilers. 

Once  program  MASTER  is  compiled,  a  run  may  be  started.  MASTER  is 
called  into  memory  and  the  run  begins.  MASTER  calls  for  BCD  input  on 
Tape  2,  and  this  is  normally  transferred  to  the  card  reader.  BCD  output 
is  given  on  Tapes  3  &  4.  Tape  4  is  normally  transferred  to  the  typewriter 
for  output  (see  BOX  10 ,  APPENDIX  V) .  The  BCD  output  on  Tape  3  is 
temporary  and  may  be  eliminated  when  the  program  is  working  correctly. 
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This  will  also  make  more  memory  space  available . 

Program  MASTER  writes  binary  output  on  Tapes  6,7,  and  8.  The 
MASTER  output  on  Tape  6  is  the  binary  input  to  program  RON .  The  MASTER 
output  on  Tape  7  is  the  regeneration  output  and  is  used  only  by  the  MASTER 
program.  The  MASTER  output  on  Tape  8  is  the  binary  input  to  program 
SLAVE.  MASTER  will  automatically  rewind  each  of  the  three  binary  tapes 
when  output  for  a  particular  tape  is  completed. 

When  MASTER  has  completed  computation  program  SLAVE  should  be 
called  from  our  system  library.  SLAVE  calls  for  BCD  input  (comment 
cards,  etc.)  on  Tape  2,  and  this  is  normally  transferred  to  the  card 
reader.  SLAVE  will  then  read  Tape  8,  compute,  and  write  its  BCD  out¬ 
put  on  the  designated  tape.  Program  SLAVE,  at  this  time,  has  PRINT 
statements  for  all  output  except  that  mentioned  in  BOX  1  of  the  SLAVE 
part  of  APPENDIX  V,  and  therefore  the  output  will  be  on  Tape  4  unless 
otherwise  directed  by  the  operator. 

When  SLAVE  has  completed  its  run,  program  RON  should  be  called 
from  the  system  library.  The  input  and  output  designations  for  RON  are 
the  same  as  those  for  SLAVE  except  that  RON  reads  Tape  6  rather  than 
Tape  8 . 

If,  for  some  reason,  MASTER  does  not  complete  its  run,  and  Tape  6 
is  not  rewound ,  program  RON  must  be  run  with  Jump  Switch  1  in  the  set 
position.  This  will  insure  that  the  correct  END  designation  (in  this  case 
a  -1)  is  placed  at  the  end  of  the  information  on  Tape  6  and  that  the  tape 
is  rewound. 
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8 .  Parameter  Adjustment 


The  two  input  parameters,  ETHRESH  and  TURN  primarily  govern  the 
cohesiveness  or  solidity  of  the  lattice.  At  ETHRESH  =  3  ev. ,  for 
instance,  some  400-500  replacements  occur  with  a  100  ev.  shot.  With 
ETHRESH  raised  to  15-25  ev.  the  number  of  replacements  drops  to  con¬ 
siderably  less  than  ten,  depending  on  the  location  and  direction  of  the 
initial  knock-on  atom .  The  best  working  values  seem  to  lie  in  the 
15-25  ev.  range,  as  preliminary  results  in  this  range  show  good  correla¬ 
tion  with  the  work  of  others .  More  accurate  values  of  ETHRESH  and 
TURN  should  be  obtained  in  the  future  from  concurrent  N-body  calcula¬ 
tions  for  use  in  the  model. 

Adjustment  of  the  TURN  parameter  also  governs  the  length  of  energy 
chains.  If  TURN  =0.0,  the  energy  chains  are  almost  non-existent. 
Atoms  vibrating  on  their  sites  are  then  re-directed  toward  their  sites 
after  their  first  collision.  Since  this  first  collision  is  usually  a  slight, 
glancing  one ,  very  little  energy  is  transferred  and  the  energy  chain  is 
not  seen.  If,  however,  the  TURN  parameter  is  set  too  high,  displace¬ 
ments  occur  among  atoms  with  energy <  ETHRESH.  Investigations  thus 
far  seem  to  indicate  that  a  value  of  about  0.1  is  the  optimum  value  for 
TURN. 

THERMAL  has  little  effect  on  the  larger  aspects  of  radiation  damage 
i.e.  chains,  channels,  displacements,  ring  displacement,  etc.,  but 
has  a  large  effect  on  machine  running  time.  By  changing  the  value  of 
THERMAL  from  0.1  ev.  to  0.5  ev. ,  the  machine  time  can  be  reduced  by 
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75%  or  more  .  The  very  low  (below  THERMAL)  energy  vibrations  are  no 
longer  seen,  to  be  sure,  but  since  this  is  not  one  of  the  areas  of  pri¬ 
mary  interest,  the  advantage  of  drastically  decreased  machine  time  far 
outweighs  the  inherent  disadvantages  „ 

When  the  penetration  of  higher  energy  (above  a  few  hundred  ev.) 
initial  knock-ons  is  the  primary  effect  under  consideration,  THERMAL 
may  be  set  as  high  as  3-5  ev.  For  initial  knock-on  energies  of  100  ev 
and  below,  THERMAL  should  be  0.1  -  0.5  ev.  to  achieve  a  reasonable 
degree  of  thermalization. 

Investigations  have  shown  that  the  radius  of  effect  (sphere  of  in¬ 
fluence) ,  SPHI,  should  lie  within  the  range  2. 5-2. 6  Angstroms  (for 
copper).  Values  much  outside  this  range  (2.2  say)  give  rather  meaning 
less  results.  We  feel,  as  a  result  of  the  studies  we  have  made,  that 

the  best  value  is  probably  slightly  less  than  the  nearest  neighbor  dis- 

l/2 

tance.  For  copper,  the  nearest  neighbor  distance  is  (2)  7  (1.807). 

We  have  used  the  value  (1.414)  (1.807)  when  not  specifically  investi¬ 
gating  SPHI  effects .  We  feel  that  any  value  larger  than  the  nearest 
neighbor  distance  will  produce  invalid  results. 

Studies  of  the  optimum  values  for  PERCENT  and  CUTOFF  are  not 
completed  at  this  time.  Preliminary  investigations  show  that  values 
above  0.1  for  PERCENT  and  below  0.1  for  CUTOFF  prevent  significant 
collisions  from  taking  place. 

The  value  of  TFAC  should  be  kept  below  1.0,  unless  the  normal 
value  of  1.0  is  used.  A  more  accurate  tracing  of  individual  particle 
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movement  is  possible  with  smaller  TFAC's,  but  machine  time  increases 


proportionately.  A  compromise  must  then  be  made  between  the  running 
time  available  and  the  individual  tracing  accuracy  desired . 

FMASS  may  be  adjusted  for  the  particular  metal  being  used.  At 
present  the  model  is  restricted  to  metals  with  face  centered  cubic  struc¬ 
ture .  The  Born-Mayer  potentials  used  are  restricted  to  copper  only 
since  values  of  the  interaction  input  parameters  (the  constants  A  &  B 
mentioned  above)  are  not  available  for  other  materials . 
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9.  Preliminary  Results 


Preliminary  results  to  date  show  that  the  existence  of  channels, 
energy  and  displacement  chains,  vacancies,  interstitials,  and  replace¬ 
ments  depend  quite  heavily,  and  quite  critically,  on  the  values  of  the 
input  parameters  used.  Before  very  satisfactory  results  can  be  obtained, 

more  reliable  values  of  these  parameters  must  be  determined. 

1  7 

A  run  corresponding  to  Gibson,  et  al's.  Run  No.  12  has  been  made 
with  a  slightly  modified  version  of  the  program,  using  a  40  ev.  primary 
knock-on.  The  primary  knock-on  was  shot  directly  at  one  of  the  atoms 
on  the  front  face  (Gibson  et  al.  start  their  primary  knock-on  inside  the 
crystallite),  at  an  angle  of  15°  away  from  the  normal  (100)  direction. 

The  values  of  the  various  parameters  for  the  run  were: 

1)  IX  x  IY  x  IZ  =  11  x  15  x  8 

2)  ETHRESH  =  20.0  ev. 

3)  THERMAL  =  0.5  ev. 

4)  CUTOFF  =  PERCENT  =  .15  (i.e.  15%) 

5)  SPHI  =2.5 

6)  TFAC  =1.0 

7)  TURN  =  0.1 

Our  results  (see  Figure  1)  are  similar  to  those  obtained  by  Gibson, 
et  al.  in  that  the  energy  chains  play  a  dominant  role  in  the  dissipation 
of  energy.  The  (100)  chain  of  our  run  is  not  as  prominent  as  Gibson, 
et  al's. ,  but  strong  focusing  is  observed  in  the  (110)  directions.  Our 
(100)  chain  was  terminated  prematurely  by  the  velocity  turn-around 
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Figure  1 .  Results  of  Preliminary  Run 
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mechanism  employed  for  lattice  cohesiveness . 
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10.  Recommendations 


Many  ideas  and  plans  occurred  to  us  during  the  development  of  the 
model  to  its  present  level  that  we  have  not  been  able  to  institute  be¬ 
cause  of  time  limitations .  These  unrealized  plans  and  ideas  are  offered 
as  possible  guidelines  for  future  development. 

We  are  convinced  at  this  time  that  the  programming  errors  in  the 
model  are  concerned  almost  wholly  with  the  LB  information  array.  As  we 
did  not  anticipate  the  great  dependence  upon  the  LB  array  when  we  began 
the  problem,  we  feel  that  the  LB  aspects  of  the  entire  model  (all  three 
programs)  must  be  overhauled,  and  more  orderly  storage  of  the  informa¬ 
tion  should  be  investigated . 

Once  this  overhaul  and  possible  re-arrangement  of  the  LB  array  has 
been  accomplished  and  the  model  is  working  as  designed,  a  more  com¬ 
plete  study  should  be  made  of  the  effects  of  parameter  adjustment. 

At  present  when  an  atom  on  the  LATER  list  doesn't  find  a  vacancy, 
a  prospective  interstitial  partner  is  chosen  from  the  MOON  list.  A 
search  is  then  made  for  vacancies  near  this  partner.  If  none  are  found, 
the  interstitial  pair  is  formed.  We  strongly  recommend  that  when  the 
prospective  partner  finds  no  vacancies,  the  rest  of  the  atoms  occupying 
nearest  neighbor  sites  should  also  be  inspected  for  vacancies  near  them. 
If,  again,,  no  vacancies  are  found  then  the  interstitial  pair  should  be 
formed  between  atom  MJ  and  its  original  prospective  partner. 

To  facilitate  this ,  we  recommend  that  the  INVAC  subroutine  be 
split  into  two  or  more  subroutines .  There  should  be  a  separate 


43 


subroutine  whose  only  function  is  to  find  the  nearest  neighbor  sites, 
not  necessarily  by  the  present  system. 

In  the  present  model,  if  a  vacancy  occurs  near  an  interstitial  pair 
after  the  interstitial  pair  has  been  formed ,  this  Frenkel  pair  will  not  re¬ 
combine  unless  the  interstitial  pair  is  hit  by  another  atom.  Provisions 
should  therefore  be  made  for  periodic  inspections  of  all  interstitial 
pairs  and  any  possible  vacancies  near  them.  A  procedure  very  similar 
to  that  discussed  above  for  making  interstitial  pairs  could  be  used. 

We  also  feel  that  subroutines  and  functions  should  be  incorporated 
to  a  much  greater  extent  than  at  present,  and  that  wider  use  should  be 
made  of  comment  cards  in  the  programs. 

A  study  of  the  relative  merit  of  a  space -wise  versus  a  time-wise 
selection  process  for  the  KHIT  list  has  been  suggested  above  and  should 
be  made. 

The  ability  to  use  atoms  of  two  different  masses  should  be  built 
into  the  model.  The  actual  use  of  two  different  masses  with  a  Born- 
Mayer  potential  must  be  delayed  until  appropriate  constants  for  the 
potential  are  found. 

ETHRESH  should  be  made  a  directional  dependent  variable  rather 
than  a  constant  as  at  present.  Since  a  function  able  to  represent  this 
directional  dependence  of  the  displacement  energy  is  not  available,  a 
table  of  thresholds  associated  with  direction  seems  the  most  promising 
course  of  action. 
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Two  memory  saving  Ideas  have  occurred  to  us: 

1)  The  possible  elimination  of  the  TLAST  concept,  and 

2)  The  possible  use  of  a  method,  already  developed,  which 
packs  the  co-ordinates  and  velocity  components  of  an 
atom  into  three  memory  cells  rather  than  the  six  required 
at  present.  This  sub-program,  SQUEEZE,  developed  by 
one  of  the  authors  (RLK) ,  saves  space  at  the  expense  of 
running  time  . 

The  use  of  these  two  ideas  would  save  considerable  memory  space,  and 
enable  the  users  to  approximately  double  the  size  of  the  available 
lattice . 

The  energy  and  velocity  scaling  processes  used  after  the  interaction 
(see  BOX  20,  APPENDICES  III,  IV,  &  V)  are  a  first  approximation  to  the 
problem  and  should  be  considered  in  that  light.  We  feel  that  at  present 
this  is  one  of  the  weakest  areas  of  the  model.  Considerable  effort  should 
be  devoted  to  this  section  in  the  future. 
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11.  Conclusions 


The  model  is  not  fully  developed  at  this  time,  but  we  feel  that 
some  valid  conclusions  may  be  drawn  from  the  results  obtained  thus  far. 
The  MASTER  program  (see  APPENDIX  III),  as  it  exists  now,  still  contains 
some  minor  programming  errors  but  we  feel  that  the  logic  is  basically 
sound. 

The  outstanding  attribute  of  the  model  at  this  time  seems  to  be  the 
greatly  decreased  machine  time  necessary  to  achieve  results  compared 
to  the  machine  time  necessary  for  an  N  -body  model.  This  decreased 
machine  time  makes  it  possible  to  use  a  lattice  as  large  as  the  machine 
memory  can  accommodate . 

We  do  not  feel  that  this  model  is  sufficiently  sophisticated  at  this 
time  to  state  definitely  the  degree  of  confidence  that  should  be  placed 
in  it.  We  do  feel,  however,  that  the  preliminary  results  show  that  with 
more  work  the  model  can  be  made  into  a  reasonable  approximation  to  the 
radiation  damage  problem.  We  anticipate  that  our  specific  recommenda¬ 
tions,  if  carried  out,  will  lead  to  a  program  which  approaches  this 
required  degree  of  sophistication. 
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APPENDIX  I 


SPECIFICATION  OF  UNITS  AND  ANNOTATED  LIST  OF  VARIABLES 
The  units  used  in  the  model  for  energy,  mass,  length,  time,  and 
velocity  are  given  below: 


ENERGY 

The  electron  volt  (ev.),  equal  to  1.60207  x  10  ^joules. 

MASS 

The  Atomic  Mass  Unit  (AMU) ,  equal  to  1 . 65983  x  10“^ 
kilograms . 

LENGTH 

The  Angstrom  unit,  equal  to  10*^  meters. 

TIME 

The  jiffy,  defined  as  10~^  seconds. 

VELOCITY 

Derived.  One  Angstrom/jiffy  =  10^  meters/second. 

The  variables  given  below  are  essentially  in  the  order  of  their 
appearance  in  the  program  listing.  However,  certain  basic  variables 
are  given  at  the  beginning  of  the  list.  An  asterisk  (*)  denotes  an  input 
variable . 


VARIABLE 

USAGE 

B  (I  ,L) 

An  array  containing  the  position  co-ordinates ,  velocity 

components,  and  energies  of  the  atoms  in  the  lattice. 

L  is  the  number  of  the  atom .  For  I  there  are  seven  cases 

1  The  X  co-ordinate  in  Angstroms 

2  The  Y  co-ordinate  in  Angstroms 

3  The  Z  co-ordinate  in  Angstroms 

4  The  X  velocity  component  in  Angstroms/jiffy 

5  The  Y  velocity  component  in  Angstroms/jiffy 

6  The  Z  velocity  component  in  Angstroms/jiffy 

7  The  energy  in  electron  volts 
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VARIABLE 


USAGE 


LB  (L) 


An  array  used  for  information  storage ,  supplementary  to 
the  B  array.  L  is  the  number  of  the  atom.  Stored  from 
right  to  left  we  find  the  following  information: 


OCTAL 

POSITION 

1  Bit  1 


Bit  2 

Bit  3 


2-5 


6-9 


10  -  13 


14 


MEANING 

1  if  the  atom  has  completed  the  current  time 
step. 

0  if  the  atom  has  not  completed  the  current 
time  step. 

1  if  the  atom  has  an  energy  >  THERMAL. 

0  if  the  atom  has  an  energy  <.  THERMAL. 

1  if  the  site  L  is  occupied  by  any  atom. 

0  if  the  site  L  is  unoccupied. 

The  number  of  the  time  step  in  which  atom  L 
exited  the  lattice.  Zero  if  L  is  in  the  lattice. 
The  number  of  the  other  atom  in  an  interstitial 
pair  with  atom  L.  Zero  if  L  is  not  a  member 
of  an  interstitial  pair. 

The  site  occupied  by  this  atom,  even  if  it  is 
the  original  site.  If  this  number  is  zero,  then 
the  atom  has  an  energy  >  ETHRESH  and  is 
wandering  through  the  lattice. 

The  number  of  the  face  through  which  atom  L 
exited  the  lattice,  called  NFACE. 
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VARIABLE 


OCTAL  (L) 
TLAST  (L) 


CUTOFF* 


PERCENT* 


USAGE 

OCTAL 

POSITION  MEANING 

15  Bit  1  1  if  the  atom  is  on  the  LATER  list  but  has  not 
been  through  INVAC . 

0  if  the  atom  is  not  on  the  LATER  list. 

Bit  2  1  if  the  atom  has  been  through  INVAC  and  has 
not  had  another  collision. 

0  if  the  atom  has  never  been  through  INVAC , 
or  has  had  a  collision  since  going  through 
INVAC. 


Bit  3  Not  used 

16  Not  used.  If  Bit  3  is  used  this  results  in  a 

negative  number  and  should  be  avoided. 
Equivalent  to  LB  (L) .  Used  for  output  purposes  only. 

The  absolute  time  from  the  start  of  the  run  to  the  time  of 
atom  L's  last  collision.  Given  in  jiffys. 

A  percentage  of  a  Time  Step  Length.  All  atoms  with  a 
time  of  collision  within  CUTOFF  %  of  a  TSL  of  the  mini¬ 
mum  time  will  take  part  in  a  near  simultaneous  collision 
with  atom  L. 

A  percentage  of  a  Time  Step  Length.  The  absolute  time 
of  completion  of  an  interaction  equals  (Time  at  start) 
plus  (PERCENT)  (TSL). 
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VARIABLE 


USAGE 


TFAC* 

A  scaling  factor  for  the  Time  Step  Lengths  . 

LULL* 

The  interval,  in  time  steps,  at  which  atoms  will  be  shot 

into  the  lattice,  provided  the  proper  jump  switches  are 

set.  Does  not  affect  operation  if  the  jump  switches  are 

not  set. 

IX,  IY,  IZ* 

The  X,  Y,  and  Z  fixed  point  dimensions  of  the  lattice, 

respectively,  such  as  9  x  14  x  15.  IX  must  be  ODD. 

A* 

One  half  the  side  of  a  unit  cell.  For  copper  this  is 

1,807  Angstroms . 

SPHI* 

The  sphere  of  influence  of  an  atom.  One  atom  has  an 

SPHI  of  (1.4 14) (A),  and  the  other  atom  is  considered 

to  be  a  point.  Given  in  Angstroms. 

ETHRESH* 

The  displacement  threshold  energy  in  electron  volts. 

An  atom  with  an  energy  >  ETHRESH  may  move  through 

the  lattice.  Atoms  with  energy  below  ETHRESH  are 

forced  to  remain  on  their  sites  in  the  lattice. 

THERMAL* 

The  thermal  energy  limit.  Atoms  with  energy  above 

THERMAL  are  allowed  to  vibrate  on  their  site.  Atoms 

with  energy  below  THERMAL  are  placed  on  site? ,  or  in 

interstitial  pairs  by  the  INVAC  subroutine . 

FMASS* 

The  mass  of  a  lattice  atom  in  Atomic  Mass  Units  (AMU) . 

JPB* 

The  interval,  in  time  steps,  at  which  the  Time  Step 

Length  is  recalculated. 
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VARIABLE 

USAGE 

MNTS* 

The  maximum  number  of  time  steps  the  program  will  run. 

MNP* 

The  maximum  number  of  particles  to  be  shot  into  the 

lattice . 

FI  (1-3) 

Used  throughout  the  program.  Usually  denotes  the 

floating  point  co-ordinates  of  a  particular  site. 

AB  (1-3) 

The  physical  dimensions  of  the  lattice  in  Angstroms  in 

the  X,  Y,  and  Z  directions,  respectively. 

NLINE 

The  number  of  atoms  in  a  line  in  the  X  direction  for  the 

undisturbed  lattice . 

NPLANE 

The  number  of  atoms  in  any  undisturbed  X  -  Y  plane  of 

the  lattice. 

LMAX 

The  maximum ,  or  total ,  number  of  atoms . 

MAXLAT 

The  total  number  of  atoms  in  the  lattice ,  excluding 

those  shot  into  the  lattice,  i.e.  the  number  of  sites. 

TEMP,  TEM 

Floating  point  variables  used  for  temporary  storage  only, 

SITE  (L) 

A  short  subroutine.  Given  the  number  of  a  site  as  an 

input,  SITE  computes  the  fixed  point  co-ordinates  of  the 

site.  These  are  then  found  in  IB  (1-3). 

TTIME 

The  total  elapsed  time  from  the  start  of  the  run  to  the 

end  of  the  current  time  step,  in  jiffys. 

OENERGY 

The  total  original  energy  of  the  lattice  in  ev.  This  is 

increased  every  time  a  particle  is  shot  into  the  lattice . 
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VARIABLE 


USAGE 


ENERGY  The  energy  in  ev.  of  the  most  energetic  atom  in  the 

lattice.  Atoms  that  have  left  the  lattice  are  not  con¬ 
sidered  in  the  computation  of  ENERGY. 

TENERGY  The  total  energy  in  ev.  in  the  lattice.  This  is  computed 

every  time  step  for  comparison  with  OENERGY. 

IENERGY  The  number  of  the  most  energetic  atom  in  the  lattice . 

MISTAKE  (l  —  1 0)  An  array  of  10  error  counters.  Their  meanings  are: 

1  The  number  of  time  steps  for  which  TENERGY  varies 
more  than  1%  from  OENERGY. 

2  The  number  of  times  a  vacancy  is  filled,  or  an 
interstitial  is  formed,  near  an  edge  or  corner  of 
the  lattice. 

3  The  number  of  times  a  negative  energy  was  computed. 

4  The  number  of  times  an  atom  formed  an  interstitial 
with  itself,  i.e.  a  false  interstitial. 

5  The  number  of  triple  interstitials  formed. 

6-9  Not  used . 

10  Set  equal  to  MISTAKE  (2)  at  the  start  of  INVAC  so 
that  the  increase  in  MISTAKE  (2)  may  be  calculated. 

NA  The  index  for  the  TSLI  and  NTSC  lists  .  Always  one 

greater  than  the  number  of  times  the  TSL  has  been 
recomputed . 
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VARIABLE 

USAGE 

NEXT 

The  number  of  the  next  particle  to  be  shot  into  the 

lattice.  NEXT  is  one  before  any  atoms  are  shot  into 

the  lattice. 

MINTS 

The  minimum  number  of  time  steps,  i.e.  the  lower  limit 

of  the  time  step  DO  loop.  Always  one  unless  the  re¬ 
generation  method  is  used. 

Til 

The  square  of  SPHI. 

MOVE 

An  indicator  used  to  determine  the  exit  point  from  sub¬ 
routine  INVAC . 

NOW 

The  lattice  site  nearest  an  atom,  computed  in  INVAC, 

used  to  determine  an  atom’s  oscillatory  center,  if  the 

atom' 3  energy  is  below  ETHRESH. 

LOCATE 

An  indicator  used  to  hold  the  number  in  Index  Register  5, 

enabling  us  to  use  Index  Register  5  without  destroying 

the  contents. 

CON  (1-4) 

Constants  used  in  the  subroutine  INTER.  These  have  no 

relation  to  the  main  program. 

W  (1-4) 

Constants  used  in  subroutine  INTER.  These  have  no 

relation  to  the  main  program. 

XC  (1-6)* 

Only  two  used.  These  are  the  constant  parameters  used 

in  the  Born-Mayer  potential ,  and  can  be  changed  to  fit 

any  one  of  the  three  types  of  Born-Mayer  for  copper, 

17 

as  devised  by  Gibson,  et  al. 
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VARIABLE 

USAGE 

Cl  -  CIO 

Constants  used  In  subroutine  INTER.  Computed  outside 

the  subroutine  as  a  time  saving  measure. 

MASK1  - 
MASK8 

\ 

Masking  quantities  used  for  logical  arithmetic  with 

LB  (L)  o  The  definition  of  the  significance  of  the  octal 

positions  in  LB  (L)  and  the  masks  themselves  explain 

ALFA 

their  specific  applicability. 

10~®  A  small  number  for  comparison  to  floating  point 

variables  to  prevent  use  of  insignificantly  small  quanti¬ 
ties  . 

ALPHA* 

The  angle  between  the  original  X  and  Y  velocity  compo¬ 
nents  of  a  particle  shot  into  the  lattice,  given  in  radians. 

BETA* 

The  angle  between  the  original  X  and  Z  velocity  compo¬ 
nents  of  a  particle  shot  into  the  lattice ,  given  in  radians . 

YENTRY  (NEXT) 

The  Y  co-ordinate  at  which  a  particle  shot  into  the  lattice 

crosses  the  X  =  0  plane  .  Given  in  Angstroms  . 

ZENTRY  (NEXT) 

The  Z  co-ordinate  at  which  a  particle  shot  into  the  lattice 

crosses  the  X  =  0  plane.  Given  in  Angstroms. 

EPB 

The  energy  in  ev.  of  a  particle  shot  into  the  lattice. 

TSL 

The  Time  Step  Length  in  jiffys.  TSL  is  calculated  so  that 

the  most  energetic  atom  in  the  lattice  will  move  a  dis¬ 
tance  equal  to  (TFAC)(A)  in  one  time  step. 

TSLI  (NA) 

A  list  of  the  first  30  TSL's  that  are  calculated.  Used 

for  output  purposes  only. 
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VARIABLE 

USAGE 

NTSC  (NA) 

A  list,  corresponding  to  the  TSLJ  list,  that  gives  the 

time  steps  in  which  the  values  given  in  the  TSLI  table 

were  calculated. 

EMAX 

A  list  of  the  values  of  ENERGY  for  the  first  30  time 

steps.  Used  for  output  purposes  only. 

N 

The  sixth  Index  Register,  and  the  number  of  the  current 

time  step.  The  current  time  step  can  always  be  found 

by  stepping  the  computer  console  and  recording  the 

value  of  this  Index  Register.  N  is  used  for  nothing  else. 

L 

The  number  of  the  atom  sequentially  under  consideration. 

The  index  of  the  DO  loop  starting  at  Statement  160. 

M 

The  number  of  the  atom  currently  under  primary  considera¬ 
tion,  usually  equal  to  L. 

J 

Usually  the  number  of  the  atom  considered  in  relation  to  M 

LATER  (NZ) 

A  list  of  atoms  whose  energies  are  less  than  THERMAL 

and  are  to  be  put  through  the  INVAC  subroutine . 

NZ 

The  index  for  the  LATER  list,  used  for  no  other  purpose. 

INMIN 

The  index  for  the  LAST  list,  used  for  no  other  purpose. 

INMAX 

The  maximum  number  of  atoms  listed  on  the  LAST  list, 

equal  to  NMAX  when  the  list  is  made. 

INPUT 

An  indicator.  If  zero,  then  the  KHIT  list  is  copied  onto 

the  LAST  list.  If  not  zero,  the  LAST  list  is  left  untouched. 
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VARIABLE 

USAGE 

NUM  (  ) 

A  list  of  all  the  atoms  in  a  box  of  side  length  3.02  (A) 

centered  at  atom  M .  Only  atoms  on  this  list  are  con¬ 
sidered  as  possible  targets  for  M.  M  is  first  on  the 

list. 

KMAX 

The  total  number  of  atoms  listed  on  the  NUM  list. 

KHIT  (  ) 

A  list  of  atoms,  derived  from  the  NUM  list,  that  atom  M 

will  hit.  Atom  M  is  first  on  the  list. 

NMAX 

The  total  number  of  atoms  listed  on  the  KHIT  list. 

TMIN  (  ) 

The  time,  in  jiffys,  of  closest  approach  between  atom  M 

and  the  corresponding  atom  on  the  NUM  list.  This 

assumes  no  deviation  from  straight  paths.  Measured 

with  the  start  of  the  current  time  step  as  a  zero  reference 

DSTANCE  (  ) 

The  distance,  in  Angstroms,  between  atorrts  when  at  their 

point  of  closest  approach  as  determined  by  their  TMIN. 

T  (  ) 

The  time,  in  Jiffys,  at  which  the  two  atoms  are  SPHI 

distance  apart.  Will  be  less  than  the  corresponding 

TMIN  if  there  is  to  be  a  collision  between  M  and  the 

atom  on  the  corresponding  NUM  list. 

ITEMP,  JTEMP ,  Fixed  point  variables  used  throughout  the  program  for 
MTEMP,  NTEMP 

temporary  storage  purposes.  No  particular  significance. 


NPAIR 

An  indicator  for  a  computed  GO  TO  statement.  NPAIR  is 

3  if  one  of  the  atoms  being  hit  is  a  member  of  an  inter¬ 
stitial  pair. 
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VARIABLE 

USAGE 

IN3 

An  indicator  that  shows  the  number  of  those  on  the  NUM 

list  that  will  be  hit  by  atom  M . 

LAST  (  ) 

A  duplicate  of  the  first  KHIT  list  comprised  for  M  when 

M  =  L.  All  atoms  on  the  list  will  become  M's  when  the 

preceeding  M  is  no  longer  able  to  have  collisions  in  this 

time  step. 

T1  -  T10 

Temporary  variables  (see  BOX  17,  APPENDIX  IV)  used  to 

speed  up  computation  time  . 

NFACE 

The  face  of  exit  of  an  atom  when  it  leaves  the  lattice. 

Faces  1,3,5  are  the  negative  X,  Y,  and  Z  faces,  res¬ 
pectively,  while  faces  2,4,6  are  the  positive  X,  Y, 

and  Z  faces  respectively. 

JA 

The  index  for  the  KHIT  list  that  shows  which  atom 

atom  M  is  colliding  with  at  the  present  time.  JA's  range 

is  from  2  to  NMAX.  Has  other  uses  in  INVAC . 

KMIN 

The  number  on  the  T  list  that  is  smallest.  T  (KMIN)  is 

the  minimum  time,  and  is  the  collision  time  for  atom 

NUM  (KMIN). 

NHIT 

The  number  of  atoms  on  the  KHIT  list  that  have  been 

ELOST 

hit  by  M. 

The  energy  in  ev.  that  M  lost  in  its  collisions  with  the 

other  atoms  on  the  KHIT  list. 
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VARIABLE 


USAGE 


SAVE  (l,K) 


COLD  (  ) 


PAR 


TIME 


DEV  (1,1-2) 


The  pre-collision  velocity  components  and  energies  of 
the  atoms  on  the  KHIT  list.  SAVE  (1,1)  is  the  initial 
X  velocity  of  KHIT  (1) .  SAVE  (2,1)  is  the  initial  Y 
velocity  of  KHIT  (1) ,  etc.  SAVE  (4,1)  is  the  energy  of 
KHIT  (1). 

A  floating  point  variable  used  in  the  energy  and  velocity 
scaling  process.  COLD  (1-6)  correspond  to  the  first 
six  co-ordinates  in  the  B  array.  They  are  the  sums  of 
the  individual  components  of  the  DEV  array  after  each 
collision,  i.e.  COLD  (1)  is  the  sum  of  all  the  DEV 
(1,1) 's  for  a  series  of  collisions .  COLD  (1)  is  then 
the  sum  of  all  the  X  co-ordinates  M  would  have  had  as 
a  result  of  the  collisions  with  the  rest  of  the  KHIT  list. 
The  scattering  parameters  in  meters ,  and  the  third  input 
to  the  INTER  subroutine. 

The  time  T  for  M's  previous  collision  with  one  of  the 
others  on  the  KHIT  list.  For  instance,  if  JA  is  3,  we 
are  considering  KHIT  (3)  in  relation  to  M„  TIME  is  then 
the  T  corresponding  to  KHIT  (2) . 

An  array,  used  in  subroutine  INTER,  holding  the  position 
co-ordinates,  velocities,  and  energies  of  the  two  atoms 
at  the  end  of  the  interaction.  The  first  index,  I,  has 
the  same  meaning  as  the  first  index  of  the  B  array. 
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VARIABLE 


USAGE 


The  second  Index  is  one  for  atom  M,  and  two  for  atom  J. 

The  B  array  is  not  changed  in  the  subroutine ,  but  the 

DEV  array  contains  the  output  of  the  subroutine.  The 

B  array  is  then  scaled  as  necessary. 

STUPID 

A  scale  factor  used  in  the  velocity  scaling  process. 

SCALE  (  ) 

A  three  member  array  used  in  the  velocity  scaling  process. 

TURN* 

An  atom  in  the  intermediate  energy  range  must  lose  TURN  % 

of  its  energy  before  it  is  re-directed  back  toward  its 

oscillatory  center. 

IT 

A  temporary  variable  used  to  facilitate  programming.  Also 

the  index  for  the  DO  loop  over  the  LATER  list. 

MANY,  NANY 

The  atoms  occupying  the  site  chosen  for  atom  M  after  M's 

energy  has  fallen  below  ETHRESH.  If  NANY  is  zero,  then 

M  and  MANY  will  become  an  interstitial  pair.  If  both  are 

non-zero,  a  "triple  interstitial"  is  formed,  while  if  both 

are  zero,  M  occupies  the  site  alone. 

JPMAX 

The  number  of  atoms  on  the  MOON  list,  i.e.  the  number  of 

nearest  neighbor  sites. 

MOON  (  ) 

A  list  of  nearest  neighbor  sites  for  a  particular  atom. 

NOW1 ,  NOW 

Temporary  variables  used  to  indicate  the  site  chosen  as 

oscillation  center  for  an  atom  in  the  intermediate  energy 

range.  NOW1  =  (NOW  shifted  27  binary  positions  to 

the  left) . 
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VARIABLE 

USAGE 

NAVY,  MAD 

Similar  to  LOCATE,  except  these  are  used  for  Index 

Registers  2  &  1,  respectively. 

ZZ  (6) 

The  time  in  seconds  necessary  for  the  interaction. 

NN,  IL,  LL 

Variables  used  to  denote  atom  numbers.  NN  is  the 

index  for  the  DO  loop  starting  at  Statement  241.  IL 

usually  denotes  a  site  number,  while  LL  usually  de¬ 
notes  the  number  of  the  "other"  atom  in  an  interstitial 

pair . 

MJ 

The  number  of  the  atom  used  as  an  input  to  INVAC . 

MM 

Usually  equals  MJ,  but  if  MJ  finds  no  vacancies  then 

MM  is  the  prospective  interstitial  partner  atom. 

I FACE  (  ) 

A  three  member  array  denoting  the  face  an  atom  entering 

INVAC  is  near  (i .  e .  within  one  lattice  unit) .  The  index 

is  1,  2,  or  3  for  the  X,  Y,  or  Z  directions,  respectively. 

NOON  (  ) 

The  list  of  vacancies  near  atom  MM  (which  is  usually  MJ) 

This  list  is  derived  by  examination  of  the  MOON  list, 

i.e.  if  an  atom  is  not  on  the  MOON  list,  it  will  never 

be  on  the  NOON  list. 

MAX 

Either  1 ,  2,  or  3.  Denotes  the  direction  (X,  Y,  or  Z  re¬ 
spectively)  of  the  axis  of  the  prospective  interstitial 

pair. 

NOTE:  This  is  the  end  of  the  list  of  variables  for  the  MASTER  program. 
Many  of  the  temporary  indices  are  not  listed  in  this  table. 
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The  symbols  listed  below  are  peculiar  to  program  SLAVE  and  the 


usage  given  applies  only  to  program  SLAVE: 


VARIABLE 

USAGE 

NTS 

The  actual  number  of  time  steps  the  run  completed. 

KET  (1-26) 

A  list  of  the  energy  ranges  for  the  histograph  distribu¬ 
tions  . 

IHISTO  (1-105) 

The  histograph  format  array. 

NBUL 

The  number  of  knock-ons  (i.e,  atoms  with  energy  above 

ETHRESH)  in  the  lattice  „ 

EBUL 

The  total  energy  of  ail  knock-ons  in  the  lattice. 

NSUBST 

The  number  of  replacements  in  the  lattice. 

NSAME 

The  number  of  atoms  on  their  original  lattice  sites. 

NVAC 

The  number  of  vacant  lattice  sites . 

INTERST 

The  number  of  interstitial  pairs  in  the  lattice. 

ETHERM 

The  total  energy  of  the  atoms  that  have  gone  through 

INVAC. 

NSIDE  (1-6) 

The  total  number  of  atoms  that  have  left  the  lattice 

through  each  of  the  six  faces  . 

ESIDE  (1-6) 

The  total  energy  of  the  atoms  that  have  left  the  lattice 

through  each  of  the  six  faces. 

ECAL 

The  same  as  TENERGY,  except  calculated  by  Program 

SLAVE. 

LMAXCAL 

The  same  as  LMAX,  except  calculated  by  Program  SLAVE 

YANGLE 

Angle  ALPHA  given  in  degrees  rather  than  in  radians. 
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VARIABLE 

USAGE 

ZANGLE 

Angle  BETA  given  in  degrees  rather  than  in  radians. 

NUMBERS  (1-50)  A  fixed  point  array  used  for  temporary  storage  only. 


IC  (1-100) 

A  fixed  point  array  used  for  temporary  storage  of  the 

histograph  channels . 

DEPTH 

The  total  penetration  of  an  atom  in  the  positive  X 

direction . 

PERDIST 

The  total  radial  penetration  of  an  atom  perpendicular 

to  the  X  direction . 

LL 

Usually  the  number  of  an  atom,  temporary  usage  only. 

XYZ  (1-3) 

Corresponds  to  the  FI  array  in  program  MASTER. 

NTOTAL 

The  sum  of  all  the  numbers  on  the  NUMBERS  list. 

NUMBIG 

The  largest  number  on  the  NUMBERS  list. 

TOTALN 

NTOTAL  in  floating  point  format. 

BIGNUM 

NUMBIG  in  floating  point  format. 

The  symbols  listed  below  are  peculiar  to  program  RON  and  the  usage 
given  applies  only  to  program  RON: 


VARIABLE 

USAGE 

T 

Temporary  floating  point  variable.  Not  to  be  confused 

with  the  T  array  in  program  MASTER. 

DX 

A  distance  parameter,  in  Angstroms.  If  any  of  an  atom's 

co-ordinates  differ  from  its  stored  values  by  DX,  the 

atom's  position  will  be  printed,  and  the  new  position 

will  be  stored . 
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VARIABLE 

USAGE 

RA 

Reciprocal  of  A. 

P  (1-3) 

The  co-ordinates  of  the  atom  under  consideration. 

IN 

The  number  of  the  time  step  the  atom  under  considera¬ 
tion  has  co-ordinates  P  (1-3). 

X  (1-3 ,N) 

Stored  co-ordinates  of  an  atom  to  be  written  out  at  a 

later  time. 

NT  (I) 

An  array  used  in  the  output  section.  The  lower  half  of 

a  storage  cell  is  the  number  of  the  atom  whose  co¬ 
ordinates  are  X  (1-3, N) ,  and  the  upper  half  is  the  time 

step  the  atom  had  these  co-ordinates. 

NO 

An  indicator.  If  greater  than  zero,  a  heading  has  been 

written  on  the  output  tape  for  the  atom  under  considera¬ 
tion.  If  less  than  zero,  no  heading  has  been  written 

on  the  output  tape  for  the  atom  under  consideration . 

S  (1-3) 

Co-ordinates  of  an  atom  in  floating  point  form  in  A  units . 

One  A  unit  is  1.807  Angstroms  for  copper.  (See  A,  above) 
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APPENDIX  II 
LATTICE  CONSTRUCTION 


0  1  2  3  4  5  6  7  8  9  10  11  12  13  14 


FIG.  II  (100)  LATTICE  PLANE,  Z  EVEN 
>Y  IX  =  15  IY  =  14 
NLINE  =  (15+1) /2  =  8 

NPLANE  =(NLINE) (IY+1)  =  (8) (14+1)  =  120 
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x  0  1  2  3  4  5  6  7  8  9  10  11  12  13  14 

T  FIG.  Ill  (100)  LATTICE  PLANE,  Z  ODD 

I - y 

To  find  the  original  co-ordinates  of  atom  77, say: 

Z  =(77-1) /NPLANE  =  76/120  =  0  Remainder  =  76 

Y  =(76) /NLINE  =  76/8  =  9  Remainder  =  4 

X  =(4) (2)  =  8 

Y  +  Z  =  9  +  0  which  is  ODD,  therefore  X  =  X  +  1  =  9 
The  original  co-ordinates  of  atom  77  were  then  (7,9,0). 


67 


APPENDIX  III 
PROGRAM  LISTINGS 

PROGRAM  MASTER 

DIMENSION  8(7 ,1500) ,IB( 15CC)  .OCTAL (1 SOC) , FI  ( 3)  .  IP (3  ), T( 20) 
DIMENSION  TMIN120) , IF ACE (3) .MISTAKE!  1 0 ) , MOON ( 1 5 ) ,NOON( 10),W(4) 

B1SIKII8S  ,c, 

DIMENSION  ZENTRY( 10) , TSL I ( 30 > , NTSC ( 30 ) , TL AS T ( 1SCC) ,LAST(2C ) 
DIMENSION  EMAX13C) ,COLD( 7) ,SAVc(4. 10 ) , SCALE! 3) 

COMMON  6, LB, OCTAL, FI , IB,T, TMIN, I  FACE , M I  STAKE , MOCN, NOON ,W , XC , XP , Z Z 
COMMON  CON.OEV, AB, A,NPLANE,NLINE ,LMAX.MASK1 , M ASK2 , M ASK3 , M A SK4 
COMMON  I  TEMP, JTEMP , MTEMP ,NTEMP , I N2,I  N3, IL ,JA  ,  JB, JC, JK, JL, JM, JP, JR 
COMMON  JT, JW, MM, KM  IN, MAX ,TEMP,NI , JPMAX,NIF, I  X. IY. I Z , NN.M A XL  A T , JJ 
COMMON  IP,TE,TEM,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C1 1,C12, T1 ,T2 
COMMON  DSTANCE,FMASS, 13, TSL 
COMMON  MASKS. MASK6,MASK7  , MASKS,  MOVE,  NCI* 

EQUIVALENCE  (  XP  »  TM I  N ) , (LB, OCTAL) 

BOX  1 

WRITE  OUTPUT  TAPE  3,  995C 
9950  F0RMAK2UH  THIS  IS  PROGRAM  MASTER  /) 

DO  9980  1=1,5 

READ  INPUT  TAPE  2,9981 

9980  WRITE  OUTPUT  TAPE  3,9981 

9981  FORMAT ( 80H 

1  ) 

BOX  2 

READ  INPUT  TAPE  2 , 803 , CUTOFF , PERCENT , TF AC , TURN 
803  FORMAT! 4F 1 0, 6 ) 

RE AO  INPUT  TAPE  2 , 1 00 1 , LLLL , I  X , I Y, IZ , A 
1001  F0RMAT(4I4,F10.6) 

READ  INPUT  TAPE  2 . 1 02 , SPH I , E THRE SH ,T HERM AL , FM ASS , JPR, MNT S , MNP 
102  F0RMAT(4F10.6,3I5) 

BOX  3 


FI  (  1)  =  IX 
F I ( 2 )  =  IY 
F I ( 3 )  =  IZ 
DO  1000  1*1,3 
1000  AB( I  )  =A*F I ( I ) 


BOX  4 

LOA( IX) , INA( 1 ) ,ENQ(0) , ARS( 1 ) ,STA(NLI  NE ),LDA( IY) . INA(  1) , 
STA  (  NPLANE  )  ,  LOA  (  I  Z )  ,  I NA  (  1 )  ,-MUI  (NPLANE)  ,STA(  LMAX)  ,STA(  MA 


, MU  I ( N L I N E  ) 
STA( MAXLAT  ) 


BOX  5 

DO  1020  L= 1 • LMAX 
CALL  SITE(L) 

1010  DO  1 C 1 5  1*1,3 
TEMP* IB ( I ) 

1015  B(  I  , L ) *A* TEMP 

ENA4(0B) ,ALS(27) ,A0D(48)  , ST A4 ( LB ) , LD A ( - 1 0. 0 ) , ST A4 ( TL AST ) 

00  1C20  1*4,7 

1020  B ( I , L )  =  0.0 

BOX  6 

WRITE  TAPE  6, I  X , I Y , I Z , AD . MAXL AT . A , SP F I , FMAS S , ETHRESH , THE RM AL , MNP 
WRITE  TAPE  6 , ( ( R ( I , J )  1  =  1,3)  J* 1  , M AX L AT ) , NL I NE , NPL ANE 
LOA(C.O) , STA (TTI ME), STA(CENERGY) , STA (ENERGY) , ST  A ( Z Z* 4 ) , EN A ( C ) 

S  T  A ( M I  STAKE* 1 ) , STA ( M I STAKE  +  2 ) , ST A ( MI  STAKE  +  3 ) . ENA ( 1 ) ,STA(  NEXT) 
STA(NA), STA (MINTS) ,LDA( SPHI ) ,FMU(SPH I ) ,STA( T 1 1 ) 

ENA ( 0 ) • STA ( M I STAKE  +  4 ) , STA ( M I  ST  AKE  +  5)  , ST A ( MOVE ) , ST A ( NOW ) 

STA ( LOCATE ) 
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ecx  7 


C  START  SET  UP  FOR  INTERACTION  SUBROUTINE 
CONI  1  )  =4. 820 78 1 099625 E-0 2 
CON ( 2  )  =  1 .089062568935F-01 
CON ( 3)=4.488872968795E-0 1 
CONI  4) *8. 6595709 1  1259E-01 

LDAI W  +  2) ,STA(W+2),LDA(W+1  ),STA(W+4) 

CC  READ  IN  POTENTIAL  PARAMETERS 

READ  INPUT  TAPE  2,  SO  *  C  XC (  I  I  1*1.6) 

50  F0RMAT(10X,6F10.5/) 

C1*l .0E-10«SPHI 

C2  =  XC I  1 )*1.602lE-l9*EXPF (XC(2)  ) 

C3*XC I  2 ) /C 1 

C4=  1 . 6598 E— 27* F NASS*  0 .25 
C5= 1 . OE-1 4 /C  1 
C6=0 . 5E 1 0*C  1 
C  7=C 1 *C 1 
C8= 1.0/Cl 

C9  =  XC 1 1  ) *1 .6021 E-19 

Cl  0=1 .0/1 FMASS*0. 51 80076  15CC). 

C  ENO  SET  UP  FOR  INTERACTION  SUBROUTINE 

BOXES. 8  AND  9 

MASK  1=7777000000 
MASK2=77770000C00006 
MASK3=77777777C0000B 
MA SK 4=7000000 00 77770B 
MASK5=1 

MASK6=100000000000000B 
MASK7=2*MASK6 
MASK8=7677777 7777 77 7770 
ALFA=1.0E-6 

USE  OF  SENSE  SWITCHES  SELECTIVE  JUMPS) 

UP, UP, UP  READ  BINARY  INPUT  OF  INTERRUPTED  PROGRAM 
UP, UP, DOWN  WRITE  BINARY  CUMP  ON  TAPE  7  FOR  REGENERATION 
UP, DOWN, DOWN  READ  IN  NEXT  PARTICLE  IF  THE  TIKE  IS  RIGHT 
DOWN, UP, DOWN  JUMP  TO  8001  WRITE  OUTPUT  AND  RETURN  TO  116 

DOWN, DOWN, UP  AT  116,  I F  SAT  I SF I  ED , GO  ON  TO  ICO,  OTHERWISE  JUMP  TO  E  ND , E  N 
UP. DOWN, UP  WILL  WRITE  TAPE  8  8 1  NARY,  BUT  NOT  THE  BCD  PART  OF  OUTPUT. 
WHILE  PRINTING  OUTPUT,  IF  WE  PUT  STOP  SWITCH  1  UP, IT  WILL  STOP  AT  THE  EN 
OF  THE  OUTPUT  TO  ENABLE  US  10  RESET  THE  PROPER  SW I TCHES, ETC . 

BOX  10 

SL  J  1  (  80  ) ,  SL  J  (  9  1  ) 

80  SLJ2I81 ) , SL J (91) 

81  SLJ3I 90) , SL J ( 91  ) 

90  READ  TAPE  7 , M I  NTS , MNT S , L M AX , B , LB , TLA 
1TSLI ,NTSC,TT I  ME, MI  STAKE, YENTRY.ZENTR 
2SPHI , FMA  SS»JPB»MNP,LAST,tMAX,VMCK»CO 
REWIND  7 

DO  9000  N=MI NTS, MNTS 

WRITE  OUTPUT  TAPE  3,8 100 ,TENERGY,N,I 
WRITE  OUTPUT  TAPE  4 , 8 1 00 , TE NERGY , N ,1 
FORMAT! FI 5. 9, 2 1 8, FI  5. 9) 

I  ENERGY=0 

THIS  IS  THE  START  OF  THE  MAIN  DO  LOO 
STEPS.  WE  CAN  CHECK  THE  CURRENT  TIM 
NOTING  THE  CONTENTS  OF  INDEX  REGISTE 
SL J 1 (82), SLJ (93) 

SL J  2 ( 83 ) , SL J ( 93 ) 

SLJ  3(93) 

WRITE  TAPE  7,N,MNTS,LMAX,B,LB,TLAST, 

1 TSL  I  .NTSC.TT I  ME, MI  STAKE, YENTRY.ZENTR 
2A,SPHI,FMASS, JPB, MNP.LAST, EMAX, VMCM, 

REWIND  7 
93  CONTINUE 


ST, TENERGY.C ENERGY, E NERGY, TSL, 
Y , NEXT , NA , EPB, LULL, I  X, IY, IZ, A, 
LC, EFOUND 


91 


8100 


82 

83 

92 


ENERGY, ENERGY 
ENERGY, ENERGY 


CF  THE  PROGRAM,  THE  ONE  ON  _ 
STEP  BY  STEPPING  THE  1604  AND 

6. 


TIME 


TENERGY.CENERGY, ENERGY, TSL, 
Y,  NEXT, N A, EPB, LULL, IX, IY, IZ, 
COLO, EFOUND 
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BOX  1  1 


SLJ  1  (  100)  ,SLJ2(85) ,SLJ( ICO) 

85  SLJ3( 100) ,SLJ(80C1 ) 

116  SLJ1 ( 9999 ),SLJ2( 9999) »  SL J3 ( ICO) ,  SLJ( 9999) 

100  LDA (NEXT), INA(-1),AJP(12C ),SLJ 1(84), SLJ ( 110) 

]?0  !  !/iji>(800l  )  ,A  JP2(  1  30  )  ,  SL  J  (  800  1  ) 

115  EN A6 ( 1 8 ) • DV l (LULL),0JP1(  130) , LDA ( NEX T ) , SUB ( MNP ) , A J P(  1 20) ,A JP2 ( 13C ) 

120  RAO(LMAX) ,LIL1 ( LMAX  )  ,  ENA  (  2  )  ,  ST  A  1  (  LB)  , LDA (-2 . 0 ) , STA 1 ( TLAS T ) 

900  REAO  INPUT  TAPE  2  »  500  , ( 8 ( I »  LMAX )  I =1  , 3 ) , B ( 7 , LM AX ) , ALP  HA, BET A 
500  FORMAT( 6F 10.6 ) 

YENTRY(NEXT)=B(2,LMAX) 

ZENTRY(NEXT)=B(3,LMAX) 

RAO( NEXT ) 

ZZ( 1 )=TANF( ALPHA ) 

FMU ( ZZ+ 1 )  STA ( ZZ  +2 ) 

Z  Z ( 3  I  =  T  AN  F (  BETA) 

FMU ( ZZ+3 )  FAC ( ZZ+2 )  FAC(l.O)  STA(ZZ+2) 
B(4,LMAX)=SGRTF(B(7,LMAX)*C10/ZZ (2)) 

B(5,LMAX)=B(4,LMAX)*ZZ(  1  ) 

B(6»LMAX)=B(4,LMAX)*ZZ(3) 

B(3,LMAX)=B( 3,LMAX)-B(6, LMAX ) /B  (  4 ,  LM  AX ) 

B(2,LMAX)=B(2,LMAX)-B(5, LMAX)/R(4,LMAX) 

125  0ENERGY=0ENERGY+B(7»LMAX) 

TENERGY=TENERGY+B(7»LMAX) 

EPB  =  B ( 7, LMAX  ) 

TSL=A*0.71  97274  588/SORTF ( ABSF ( B ( 7 . LM AX ) /FMASS ) ) *TF AC 
+LDA ( NA ) , INA(-30) , +  AJP ( L  +  1 ) , A JP2 ( 1 40 ) , L I L 1 ( NA ) , LDA ( TSL  ) 

ST  A  1 (TSL I ) ,ENA6(0) , STA 1 (NTSC) , RAC (NA ) ,SLJ( 140) 

BOX  12 

130  ENA6( 0) , ENQ( 0 ) ,DVI ( JPB ) , CJP1 ( 1 40  ) 

135  TSL=A«0. 71 97274588/SORTF (ABSF (ENERGY /FMASS) )«TFAC 
+LD A ( NA ) , INA(— 30) , +AJP ( L+ 1 ) , AJ  P2 ( 140 ) 

LI  LI (NA) ,LDA( TSL) ,STA1  (TSL I ) , ENA6(0)  , ST  A  1 ( NTSC ) , RAO ( N A ) 

BOX  13  . 

140  LDA ( TTI ME ) »FAD(TSL)»STA( TTIME) 

LDA ( OENERGY ) ,FMU( 1 .01 ) ,FSB( T ENERGY) ,  A JP ( 155 ) , AJP3 ( 1 50 ) ,LDA ( 0 .C ) 

FSB  (OENERGY)  .FMU  (  0.99.)  ,  FAD  (  T  ENERGY ),  A  JP  (  1  55 )  ,  A  JP2  (  1  55  ) 

150  RAO ( M I STAKE+  1  ) 

BOX  14 

155  EN A6 ( 0 ) , I NA (-3 1 ) , +A JP2 ( L  +  3 ) , LD A ( ENER GY ) , I N I  6 ( - 1 ) , STA6 ( EM AX  ) 

♦  INI6( 1 ) ,+LDA (0.0) , STA (ENERGY) ,STA(T  ENERGY) 

DO  160  L=1,LMAX 
TENERGY=TENERGY+B(7»L) 
fLDA4(LB) ,SCL( MASK5)  , ST  A4 ( LB ) 

160  CONTINUE 

BOX  15 

C  THIS  IS  THE  START  OF  THE  MAJOR  DO  LOOP  CN  THE  PARTICLES 
DO  8000  L= 1 , LMAX 

WITHIN  A  PARTICULAR  TIME  STEP,  THIS  IS  THE  START  OF  THE  MAIN  CC  LOOP  OV 
ALL  THE  PARTICLES  IN  THE  LATTICE.  WE  CAN  TELL  WHICH  PARTICLE  IS  BEING 
SIDERED  BY  STEPPING  THE  16C4  AND  CHECKING  INCEX  REGISTER  4. 

LDQ4 ( LB ) , LDL ( 2B ) , A JP ( 8C0C ) , L DL ( 77771 E1.AJP1 (8000) 

LDL (MASK7)»AJP1( 8C00) 

165  DO  161  K= 1 »  50 

161  LATER( K )=C 

ENA ( 1 ) ,STA(NZ) ,ENA( 2) »STA( IN MIN) ,ENA (0) »STA( INPUT ) 

DO  164  K=1 ,20 

ENA(0),STA3( NUM ) ,STA3(KHIT),STA3(LAST) , LCA(C.O) ,STA3(T) 

STA3 ( DST  ANCE ) ,STA3(TMIN) 

164  CONTINUE 

.  S I L  4 ( ITEMP),LIL5( ITEMP) 

C  THIS  IS  THE  END  OF  THE  USE  OF  L  AS  SUCH  FOR  THE  REST  OF  THE  CC  LCOP 
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BOX  16 


200  ENA(-C) v STM  IN 3) , ENA5 ( 0 ) , A JP ( 277 ) , ENI  3 ( 2 )  ,SIL5(NUM*1 ) ,SIL5(KHIT+1) 
L0Q5 (LB) ,LDL(MASK7) . A J  P 1 (277) 

NPA  I  R= 1 
DO  201  1=1,3 

201  T(I1=B(I,M) 

LD02?LB )  *  A JP 1  ( 20  3 ) 

DO  2C2  1=1,3 

IF(ABSF(B(I»J)-TII ) > - 1 . 5  1  * A  )  202,20  3,20  3 

202  CONTINUE 

SIL2(  I  TEMP) ,LAC( I  TEMP) , I  NAS ( 0 ) ,AJP(203) ,ENA2(0) ,STA3(NUM ) ,INI3( 1  ) 

203  CONTINUE 

I NI 3 (- 1 ) , ENA3 ( 0 ) ,STA(KMAX) , INA<-2) ,A JP3(277) ,ENA<0) ,STA( IN3) 
STA(NFACE) 

BOX  17 

DO  215  K= 1 , KMAX 
LDA3 ( NUM )  ,STA( I  TEMP ) 

T 1 =  B ( 1 , I  TEMP )-B (  1  , M ) 

T2=B(2,ITEMP)-B(2,M) 

T3=B( 3, ITEMP)-B(3,M) 

T4=B ( 4 , ITEMP )-B( 4,M) 

T5=B(5, ITEMP )-B( 5  »M ) 

T6=B(6, ITEMP) -B (6, M) 

T7=T4*T4+T5*T5+T6*T6 
T8  =  T1*T1 +T2*T  2+T  3*  T3 
LDA ( T7) , FSB (ALFA) ,A JP3(2 14 ) 

T9=T1*T4+T2*T5+T3*T6 

LDA (  1.0) , FOV ( T7 ) ,  STA ( T 10 ) , L AC ( T9 ) , FM U ( T 1 0 ) , S T A 3 ( TM I N ) 
DSTANCE(K)=SQRTF(ABSF( (T1+T4*TMIN(K)  )«(T  1  +  T4*TMIN(K) )  +  (T2  +  T5*TPIN( 
IK) )«(T2+T5*TMIN(K) )+(T3+T6»TMIN(K) ) * ( T 3+ T6* TM IN ( K ) ) ) ) 

F  SB ( SPH I ) ,AJP(214) ,AJP2( 214) 

T(K)=TMIN(K)-SQRTF(ABSF(TMIN(K)*TMIN(K)-(T8-T1 1 )/T7) ) 

211  LD A3  (  T  )  , FSB  (  TSL )  ,  AJ P(  2  1 4  ) ,  A  JP2  ( 2  14  ) .  FAC  (  TSL  )  ,F AD (  TSL  )  ,  AJ P3  (  2  1 4  ) 
L0A3(T) ,FDV(TSL) ,FSe( 0.95999 ) ,AJP( 21  2) ,AJP2(214) 

212  LIL1 ( ITEMP) ,L0A3(T) , FAO ( TT I  ME ) ,FSB ( T SL ) , FSB  1 ( TLAST ) , AJP(  2 1 4 ) 
AJP3(214) ,LDA3(T) , FAD ( TT IME ) , F SB ( TSL ) , FS B5 ( T L AST ) , A JP  (  2 1  4 ) 

A JP  3 ( 2 1 4 ) 

213  RAO ( IN3 ) , SL J ( 2 1  5 ) 

214  LDA(TSL) ,FMU( 1 . 0 E+2 ) , STA 3 ( T ) 

215  CONTINUE 


BOX  18 

LDA( IN3) »AJP1(216),SLJ(2?7) 

216  LOA (KMAX) .INA(-1 ) , STA( JA  )  ,LCA( TSL) ,FMU( 1 .CE+3) , STA (TEMP) 
00  217  K=2 , KM AX 

^  LDA3 ( T ) ,FSB( TEMP) ,AJP2 (217), LDA3 ( T ) ,  STA ( TEMP ) , SI L3 ( KM  IN ) 
2fl7  CONTINUE 

L I L 1 (KMIN) ,ENI2( 2) ,LDA (TSL) ,FMU( CUTO FF ) , STA ( T1 ) 

00  219  K= 1 ,KMAX 

LDA3 ( T ) , FSB  1 ( T ) , F SB ( T 1 ) , A JP ( 2 1 8 ) , A JP 2 ( 2 1 9 ) 

218  L0A3 ( NUM ) ,STA2(KHIT) , INI21  1  ) 

219  CONTINUE 

I N I  2 (  —  1 ) , SIL2(NMAX) ,LDA( INPUT) ,AJP1(  301  ) ,SI L2( INMAX) 

DO  300  K= 1 , INMAX 
LDA3(KHIT) ,STA3(  LAST) 

300  CONTINUE 

301  RAO ( INPUT) ,LOA(NMAX) , INA (-1 ) ,AJP1 ( 22 0 ) , SL J ( 2 77 ) 

220  ENA(2) »STA(JA) , ENA ( 0 ) , ST  A ( NH IT ) , LOA( 0.0 ) , ST  A ( ELCST ) 
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BOX  19 


00  451  K=  1 , NMAX 
J=KH I T ( K ) 

DO  351  1*1,4 
351  SAVE( I,K)*B(I+3, J) 

00  43C  1=1,7 
430  COL D( I ) =0.0 
221  J=KH I T ( J  A ) 

224  L0Q2 ( LB) ,LDL ( MASK  1 ) ,+AjP( L  +  2 ) , EN  A( 3)  , -ST  A  (NPA  IR  ) 
00  245  K= 1 » KMAX 

LAC3( NUM) , INA2 ( 0 ) ,AJP(246) 

245  CONTINUE 

246  PAR=CSTANCE(K)*1 .0E-10 
DO  240  1=1,3 

B( I , J)=B( I , J) +B( 1+3, J)*T (K) 

LDA  (  J  A  )  ,  I N  A  C  -  2- )  ,  AJPI225) 

B(  I  »  M ) *B ( I , M ) +B ( I  +  3 ,M ) * ( T ( K ) -T I  ME ) 

SL J ( 240) 

225  B( I ,M)=B( I,M)+B( I+3,M)*T(K) 

240  CONTINUE 

T I ME  =  T ( K ) 

LDA2ILB) ,SCL( MASK7) ,STA2 (LB) 

380  CALL  INTERIM, J, PAR) 

TLAST( J)=TTIME-TSL+TIME 
805  DO  227  1=1,7 

COLD ( I ) =COLD ( I ) +DEV( 1,1) 

227  B ( I , J ) =DEV (1,2) 

COLD ( 7 ) =0.0 

EL0ST=B(7, J)-SAVE(4,JA)+ELCST 
SL J ( 22 1 ) 
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BOX  20 


T2=B(7,M)-ELOST 

AJP ( 356) ,AJP2(229) , LDA (0.0) , STA ( TEMP ) , S T A( T EM ) 

356  DO  357  1=4,7 

357  B(  I,M)=0.0 

LDA ( T2 ) , A JP ( 229 ) 

DO  358  K= 1 , NMAX 
J=KH  IT ( K  ) 

520  DO  352  1=1,3 
TEMP=TEMP+SAVE( I , K ) »S  AVE ( I ,K) 

TEM=TEM+B( 1+3, J) »e( 1+3, J) 

352  CONTINUE 

358  CONTINUE 
TEMP=TEMP/TEM 
STUPID=SQRTF( ABSF(TEMP) ) 

DO  355  K=2 , NMAX 

J=KH I T ( K ) 

521  DO  354  1=4,6 

354  B( I , J)=STUPIO*B( I,J) 

B ( 7 , J)=B(7,J)»TEMP 

355  CONTINUE 

229  TEM=NHI T 
TEM=1 .O/TEM 
DO  359  1=1,3 

359  B(I»M)=COLD(I)»TEM 

LDA ( T2 ) , A JP (231 ) , A JP3 ( 2 3 1 ) , LDA ( S AVE+ 4 ) , FSB ( ETHRE SH ) , A JP2 ( 503 ) 
LAC  (  SAVE  +  4)  ,  FMU  (  TURN  )  ,  FAC  (  ELOS.T )  ,  A  JP  (  L  +  2  ) ,  A  JP2  (  L+ 1  )  ,  SL  J  (  50  3  ) 
LDQ51LB) , LOL ( MASK2 ) ,ARS( 27) , STA( IL ) , LDA( 0.0 ) , STA ( TEM) 

CALL  SITE  (IL) 

DO  501  1=1,3 
TEMP=I B ( I ) 

SCALE( I )=TEMP*A-B( I .M) 

FMU 1 ( SCALE ),FAO( TEM), STA (TEM) 

501  CONTINUE 

LDA(TEM) , FSB ( ALFA), AJP(5C3), AJP3( 503) 

TEM=SCRTF(ABSF(T2*C1 O/TEM) ) 

DO  502  1=1,3 

502  B( I+3,M)=TEM*SCALE( I) 

B ( 7 , M  )  =T2 

GO  TO  231 

503  LDA(NHIT) »INA(-1 ),AJP1 (505) 

DC  504  1=4,7- 

504  8 ( I ,M)=DEV( I,  1 ) 

GO  TO  231 

505  TEMP=0 .0 

DO  329  1=4,6 

LDA 1 (COLD)fFMUI(COLD) , FAC ( TEMP ), STA ( TEMP) 

329  CONTINUE 
B ( 7  ,  M  )  =  T2 

TEMP=SQRTF(ABSF(C10»B(7, MI/TEMP) ) 

230  B( I ,M)=COLD( I )»TEMP 
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BCX  21 

231  ENA  (0),  ST  A  (NX),  STM  NY)  , +  LDA(NHIT)  ,  AJ  P  1  (L  +  2)  ,  STA(  IN3)  ,  SLJ1277) 
TLAST(M)=TTIME-TSL+TIME+TSL*PERCENT 
LOC  ATE  =  0 
TEMP=B( 7, M) 

FSB( ETHRESH) , A JP ( L+9 )  ,  +  A JP3 <  L  +  8 ) , LOO ( M ASK2 ) , LCL5 ( LB ) , AJP (L  +  4  ) 
+ARS(3),STA( I  TEMP) ,LIU1 ( ITEMP) ,LDA1(  LB) ,SCL (LB) ,STA1 (LB) ,LDA5(L8  ) 
SCL( MASK2 ) ,SST( 26) , STA5( LB) , +SLJ (514 ) .+LDA( TEMP ), FSB ( THERMAL ) 

A JP2 ( L+1 ) , SLJ ( L  +  3) ,LDG(MASK2 ) ,LDL5(LB) , A JP 1 ( 5 1 3 ) , SL J ( L+2 3 ) 
+LDA5(LB),+SCL(2B),STA5( L 8 ) , LOC ( M ASK  2 ) , S TL (  I  TEMP)  • 

LDQ ( MASK  1 ) ,STL( J TEMP) ,LOA( JTEMP) , AJP 1 (L+6),+LDA(  I  TEMP)  ,AJP( L+4  ) 
ARS ( 27 ) , STA( I  TEMP ) , LI L 1 (  I  TEMP ) ,L0A1(  LB ) , SCL ( 4 B ) , ST A  1  ( LB )  , SL J ( L  +  4  ) 
+  ARS(  15),STA(  JTEMP)  ,LIL1  (  JTEMP  ),LDA1  (  LP. )  ,  SCL  (MASK1  ),STA1(LB) 

LD A5 ( LB ) »  SCL ( MASK3 ) , LDQ ( M A SK6 ) , STL ( I T EMP ) . SC L ( MASK 6 ) , ST A  5 ( L P  ) 

LD A ( I  TEMP) , A  JP 1 ( L  +  4 ) , LI  LI ( NZ ) , ENA5(0 ) , STAl { LATER ) ,  IN  I  1  ( 1  ) 

+  SIL1  (NZ) »  +LDA ( T  EMP )  , A JP ( L  +  2 ) , AJ P2 (L  +  1  ) , R AO ( M I  ST AK E  +  3 ) , E NA ( C ) 

STA ( LOCATE ) 

600  SIL51MJ) ,ENA ( 1 ), STA ( MOVE ), ENA ( 0 ), STA ( NOW ) 

CALL  INVAC ( M J  ) 

ENA(C) , STA (MOVE) , LO A5 ( LB ) , SC L ( MA SK7) ,  ST A5 (  L P ) ,L CC ( NO W ) , Q J P ( 5  1  C  ) 
QLS ( 27 ) »  STQ ( NOW  1  ) 

509  LDA5( LB ) ,SSU(MASK2) , STA5 ( LB ) , L I L 1 ( NOW ) ,LD01  ( LB) , LDL ( 48 ) , AJP 1 ( L  +  3  ) 
-L0A1  (LB)  »SST(4B)»STA1  (LB) ,SLJ(513),  +ENA(0) , ST A ( MANY ) , ST A( NANY ) 

00  507  I  =  l ,  L  M  AX 

LOOK  LB)  ,LOL  (MASK2  )  ,  SUB  (  NOW  1  )  »  AJP  1  (5  C7  )  ,LOL(  MASKl  )  ,  AJP(  508) 

ARS ( 1 5 ) »  STA  (  MANY ), LOO (NANY) ,QJP1 (507 ), STAl NANY) 

507  CONTINUE 

508  ENA1  (0) , AJP1 ( L+4) .  +  LIL1 ( NOW )  , LCA 1 ( LB ) , SCL ( 4 B ) , ST  A 1  ( L6),LCC(NCWl ) 
SLJ15C9) ,LDA( MANY), SUP (NANY)  ,  +  AJP1 (L  +  1  ) , RAO ( M I S T AK E+ 5 ) ,LDA ( NANY  ) 
ALS  (  15) , STA ( I  TEMP) ,ENA5(C) ,ALS( 15) ,STA( JTEMP ) , LO A5 ( LB ) , LOO ( M ASK  1  ) 
SSU( I  TEMP), STA5(LB) ,  +  LILl (NANY) , +LOA 1 ( LB ) , SSU ( JT EMP ) , ST A  1 ( L B ) 

ENA  1(0), SUB (MANY) ,AJP(51 3 ) , L I L 1 ( MANY ) , SL J ( L- 3 ) 

510  DO  51 1  K  = 1 , JPMAX 

L0A31M00N) ,STA( I L) , LD A ( 0 . C ) , ST A3 ( OST ANCE ) 

CALL  SITE  ( I L  ) 

DO  511  1=1,3 
TEMP=  I B ( I  ) 

TEMP=B( I ,M)-TEMP*A 

511  DSTANCE(K)=DSTANCE (K)+ TEMP* TEMP 
1  =  1 

DO  512  K=2, JPMAX 

L0A3 (DSTANCE ) ,FSB1 ( OST ANCE ) , AJP3 ( L+1  ) , SL J ( 5 1 2 ) , SIL3 ( IT ) , L I L 1 (  I  T ) 

512  CONTINUE 

L0Q1 (MOON) ,STQ( NOW) ,0L S ( 27 ) , STC ( NCWl ) . SL J ( 509 ) 

LOA  (  LOCATE  )  ,  A  JP  1  (  5 1  6  )  ,  LD  A  (  T2  )  ,  A  JP  (  L+ 2  )  ,  A  JP2  (  L+ 1  )  ,  R  AO  (  M I S  T  AK  E  + 3  ) 
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BOX  22 


514  ENI  1  (2) 

515  LIJA1  jKHIT  )  ,STA(  IT) ,  LIL2(  IT) 

FSB? ETHRESH) , AJP (L  +  9)  , +A JP 3 ( L*8) ,LDQ (MASK2 )  .LCA2 (LB) , STL ( I  TEMP  ) 

SCL ( MASK2 ) »SST(2B) , STA2 ( LB ) , LOA ( ITEMP),AJP(233),ARS(3),LIU3IITEMP) 
LDA3 ( LB ) * SCL ( 4B ) ,  ST  A3  (  LB  )  , -*-S  L  J  (  2  3  3  ) , +LCA( T1 ) , FSB (THERMAL ) 

AJP(L+1)  ,AJP2(L+17) ,LDC(MASK2) ,LDA2( LB) , SCL ( 2fi ) , ST A2 ( LB ) 

STL (  I  TEMP) ,LDG(MA$K1  ) t  STL ( J TEMP) , LDA ( J T EMP ) , A JP ( L +8 ) , ARS (15) 

ST A ( JTEMP ) , L I L3 ( JTEMP ) ,LCA3( LB) , SCL(  MASK  1 ) , STA3 ( LB ) , LDA2 ( LB ) 
SCL(MASK3),SST(MASK6) , STA21LB) , L I L3( N Z ) , ENA 2 ( 0 ) , ST A3 ( L AT ER ) 

RAO(NZ) •  SL  J ( 2  33 ) , LDA ( I  TEMP ) , ARS ( 27 ) ,  STA (  I  TEMP ) ,L I L 3 ( I  TEMP ) 

LDA3(LB) ,SCL(4B) ,  +  SLJ (L-E) ,♦ LOG (MASK  2) , LDA2 ( LB ) , SST ( 28 ) , ST A2 ( LB ) 
LOL2 (LB) ,AJP1 (233) ,SIL5( LOCATE) ,SIL2 ( N A V Y ) , L I L5 ( NAVY ) ,SIL1 (MAO) 
SLJ( 600) 

516  LIL2(NAVY) ,LI L5( LOCATE) , LIL1 (MAD) 

233  ENA  1 (0) ,SUB(NMAX) ,A JP ( 234  ) , I NI 1 ( 1 ) ,S L J ( 5  15 ) 

BOX  23 

234  LOA ( NHIT) , AJP1 (241 ) , LDA( C.O) , STA ( ZZ+ 6 ) , SL J ( 283) 

241  00  282  NN= 1 « NMA X 

J=KH I T ( NN ) 

DO  280  K=1,KMAX 

LAC  3 ( NUM ) » INA2( C ) »A JP ( 28 1 ) 

280  CONTINUE 

281  LDA(NN) , I NA (  —  1 ) ,  +  A JP ( L+2 ) , LO A3 ( T ) , -S T A ( T I  ME ) 

DO  235  1=1,3 

235  B(  I  ,  J)=B( I , J)-( TIME  +  ZZ( 6 )*1 ,0E+14)*B ( 1+3, J) 

282  CONTINUE 

BOX  24 

283  GO  TO  (277,277,236,277) ,NPAIR 

236  DO  238  K=2,NMAX 
J=KH I T ( K ) 

T2=B(7, J) 

FSB (ETHRESH) , A JP ( 52 6) , A J P3 ( 5 26 ) , LDG2 ( LB ) , LDL ( M A SK2 ) , AJ P( 2 3 8 ) 

STA (  IL) , LDL (MA SKI ) , AJP ( 2  38 ) , STA ( LL ) ,  LDA2 ( LB )  , SC L ( MASK3 ) , STA 2 ( LB ) 
LDA(LL) ,ARS( 1 5) ,STA(LL) 

T4=B(7,LL ) 

FSB (ETHRESH) ,AJP3(L  +  7) , L I L 1 ( LL ) , LDA1  ( LE ) , SCL ( MA SK3 ) , S T A 1 ( LB ) 

LDA (  IL), ARS( 3) , S T A (  I  L  ) , L  IU 1 (  I L ) , LD A1  ( L B ) , SCL ( 48 ) , S T A  1 ( LB ) , SL J ( 23 e ) 
LIL  1  (LL)  ,LDA1  (LB)  ,SCL  (M.ASK1  )  ,STA1  (LB  )  ,  LDQ(  MASK7  )  ,  STL  (  IT)  ,LDA(  IT) 
AJP (238) ,LDA(  IL) , ARS (27)  ,  +  STA( IL) 

CALL  SITE  ( IL) 

DO  237  1=1,3 
TEMP= I B ( I ) 

237  B ( I , LL ) =A*TEMP 
GO  TO  238 

526  LDQ2 (LB) ,LDL( MASK1 ) ,AJP(  238) ,ARS(  15)  , ST A ( LL ) , L I L 1 ( LL ) 

T4=B( 7, LL ) 

FSB ( ETHRESH) , AJP 3 (238) , LCA1 ( LB ) , SCL( MASK3 ) , STA 1 (LB) ,LDA2 (LB ) 

SC  L ( MASK  1 ) ,STA2 (LE) 

238  CONTINUE 

BOX  25 

277  LDA ( IN3) ,AJP1 (200) , LDA ( I NMIN) , SUB( IN  MAX) , AJP (302 ) , AJP2(  30  3  ) 

302  M=LAST( INMIN) 

RAO ( INMIN)„LDQ5(LB) ,LDL(MASK6) ,AJPl(  277) ,SLJ(20C) 
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BOX  26 


303 


304 

525 

350 

305 


261 

262 
2  50 

251 

252 

254 

255 
307 


331 

7998 

8000 


t  309 
310 


31  1 
312 


314 

315 
257 


332 

320 

316 


9000 


DO  3C5  J  =  1 »  I NM  AX 
K  =  L  AST ( J ) 

LDQ3 ( LB ) ,  LDL ( MASK6 ) ,  AJP1 (52b ) 

DO  304  1='  3 
B ( I f  K ) =B ( 

LDAi ( LB ) f 
I F ( ENERGY 
ENERGY=B ( 

I EN  ERGY  =  K 
CONTINUE 

BOX  27 

LDA(NZ), INA(-l) »  A JP ( 262 ) ,AJP3(262),STA(IT1) 

DO  261  I  T= 1 « I T 1 
M J  =L A  T  ER ( IT) 

CALL  INVAC(MJ) 

ENA(O) ,STA( IU)*SIU6(IU)  LDA(MJ)  RAD(IU) 

WRITE  TAPE  6»IU»(B(I*MJ)  1  =  1 , 3  ) 

CONTINUE 

'  BOX  28 

DO  7998  J  =  1 , INMAX 

LDA2 (LAST) ,STA( I  TEMP) , LIL5(  l  TEMP  )  ,  LD  Q5  (  L  B  ) ,LDL( MASK7 ) , AJP1 ( 7999) 
DO  252  1=1,3 

IF ( B ( I ,M) -AR (  I  )  )  251.251,254 
I F ( B (  I , M ) )  255,252,252 
CONTINUE 
GO  TO  307 

ENA  1(0),ALS(1  )»STA(NFACE)*SLJ(252) 

ENAl  ( 0 ) , ALS (  1  )  , I NA (- 1  ),STA(NFACE)  ,SLJ(252) 

LDA( NFACE) ,AJP( 7998 ) «  ALS ( 36 )  ,  INA6(0)  ,ALS(3),STA(ITEMP),LDA5(LB) 
LDQ ( MASK4 ) ,SSU( I  TEMP ) , LDC ( MA SK2 ) , STL  ( I  TEMP)  , LDQ( MASK  1 ) ,STL ( JTEMP  ) 
SCL (MASK3 ) ,STA5 (LB) ,ENA( C) ,STA( NFACE ) , LDA( JTEMP) , AJP 1 ( 33 1 ) 

LDA( I  TEMP) ,ARS(27),STA( I  TEMP)  ,LIL1 (I  TEMP) ,LOAl (LB) ,SCL(4B) 

STA1  (LB) , SL J ( 7998 ) 

ARS(  15) ,STA( JTEMP) , L I L 1 ( JTEMP) ,LDAl  (  LB)  , SCL ( MASK  1 ) , STAl ( LB) 

CONTINUE 

CONTINUE 

THIS  IS  THE  END  OF  THE  MAJOR  DO  LOOP  ON  THE  PARTICLES 
BOX  29 

ENA ( 0 ) , ST  A ( IT ) ,SIU6( IT) 

DC  9C00  L  =  1 , LM AX 

LDQ4 (LB) ,LDL( 2B) , AJP (309 ) ,SIL4( IT) 

WRITE  TAPE  6,IT,(R(I,L)  1=1,3) 

LDQ4(LB) ,LDA(MASK8) ,STL4(LB) ,LDL(77771B) ,AJP1 (9000) 

LDL ( MASK7) ,AJPl ( 320 ) 

DO  310  1=1,3 

B(I,L)=B(I,L)+B(I+3,L)*TSL 
DO  312  1=1,3 

IF ( B( I , L ) — AB ( I ) )  311,311,314 
I F ( B (  I , L ) )  315,312,312 
CONTINUE 
GO  TO  257 

ENAl ( 0) ,ALS( 1 ) ,STA( NFACE ) ,SLJ(312) 

ENAl (0) ,ALS( 1 ) , I NA ( —  1 ) , ST A (NFACE ) , SL J ( 3 1 2 ) 

LDA ( NFACE ), AJP (320) , ALS ( 36 ) , I NA6 ( 0 ) , ALS ( 3 ) , S TA ( I TEMP ) , LD A4 ( LR ) 

LDQ ( MASK4 ) , SSU ( I  TEMP)  , LDC ( MASK 2 ) , STL ( I  TEMP)  ,LDQ(MASK1),STL( JTEMP) 
SCL ( MASK 3 ) ,STA4(LB) ,cNA(C) , STA (NFACE ) ,LDA( JTEMP ) ,AJPK  332) 

LDA(  I  TEMP ) , ARS127)  ,STA( I  TEMP) ,LIL1  (I  TEMP)  ,LCA1  (LB) ,SCL(4B  ) 

STAl  (LB)  ,  SL  J  (  320  )' 

ARS  (  15)  ,  STA  (  JTEMP)  ,LIL1  (  JTEMP.)  ,LCA1(  LB),  SCL  (  MASK  1  )  .STAl  (  LB) 

I F ( ENERGY-B ( 7 , L ) )  3 1 6 , 90C C , 9000 
ENE  RG Y=B ( 7, L ) 

I ENERGY=L 
CONTINUE 

THIS  IS  THE  ENC  OF  THE  MAJOR  DC  LOOP  ON  THE  TIME  STEPS 


-B ( 7 , K ) )  350,305, 


7 ,  K ) 


305 
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BOX  30 


F 1 3. 7,7H  JIFFYS  ) 
F13.7,7H  JIFFYS) 


8001  WRITE  TAPE  8 ,  I  X ,  I  Y , I Z , AB , N'L I NE , N PL AN E , MAXL A T , LM A X , A , S PH  I , FM A S S , 

1 E THRESH, THER MAL, EPB, LULL, ALPHA, 8 ETA,  N  ,MNTS»TTIME»TSLI»NTSC» 
20ENERGY,TENERGY, MISTAKE,  YENTRY, ZENTRY ,N A , EMAX, EFOUND 

WRITE  TAPE  8, ( LB ( L ) , (R( I ,L)  1=1,7)  L=1,LMAX) 

REWIND  8 

♦  SL  S  1  ( L  + 1 ) ,+SLJl (L  +  l ) , SL J ( L  +  2 )  ,SLJ2( L  +  1  )  ,SLJ3(999B) 

WRITE  OUTPUT  TAPE  3,  9951,TIIME 

9951  FOR  MAT ( 1 H 1 , 1 9H  THE  TOTAL  TIME  IS 
WRITE  OUTPUT  TAPE  3,  9952, TSL 

9952  FORMAT ( 25H  THE  TIME  STEP  LENGTH  IS 

I F ( N  )  9974,9974,9973 

9973  WRITE  OUTPUT  TAPE  3 ,99 53 , ENERGY , N 
GO  TO  9972 

9974  WRITE  OUTPUT  TAPE  3 , 9953 , ENE RGY , MNTS 

99530F0RMATI9H  ENERGY=  F15.7,15H  ELECTRON  VOLTS/14,  27H  TIME  STEPS  HAV 
IE  BEEN  TAKEN  ) 

9972  WRITE  OUTPUT  TAPE  3 , 9954 , T EN ERGY , OEN ERGY 

99540 FORMAT! 2 1H  THE  TOTAL  ENERGY  IS  F12.7.15H  ELECTRON  VCLTS/24H  THE  C 
1RIGINAL  ENERGY  IS  F12.7.15H  ELECTRON  VCLTS) 

I F ( M  I  STAKE ( 1 ) )  9955,9955,9956 

9955  WRITE  OUTPUT  TAPE  3,9957 

9957  FORMAT ( 30H  THE  ENERGY  CHECK  DID  NOT  FAIL  ) 

GO  TO  9965 

9956  WRITE  OUTPUT  TAPE  3 , 99 59  ,  M I  ST AK E ( 1 ) 

9959  FORMAT ( 24H  THE  ENERGY  CHECK  FAILED  I4,16H  DIFFERENT  TIMES) 

9965  IF ( MI  STAKE ( 2 ) )  9961,9961  ,9962 

9961  WRITE  OUTPUT  TAPE  3,9963 

9963  FORM  AT ( 62H  NO  INTERSTITIALS  OR  VACANCIES 
1C0RNER  ) 

GO  TO  9805 

9962  WRITE  OUTPUT  TAPE  3 , 9960  ,  M I S T AKE ( 2 ) 

9960  FORM  AT ( 60H  AN  INTERSTITIAL  CR  VACANCY  CCCURRED 
1RNER  14, 1 6H  DIFFERENT  TIMES) 

9805  I F { M I  STAKE ( 3 ) )  9800,9800,980  1 

9800  WRITE  OUTPUT  TAPE  3,9803 

9803  FORMAT ( 36H  NO  NEGATIVE  ENERGIES  WERE  COMPUTED 
GO  TO  107 

9801  WRITE  OUTPUT  TAPE-  3 , 9804  ,  M  I  S  T  AKE  (  3 ) 

9804  FORMAT ( 3 1H  A  NEGATIVE  ENERGY  WAS  COMPUTEC  14, 

1  /) 

107  I F ( M I  STAKE ( 4 ) )  170,170,  172 

170  WRITE  OUTPUT  TAPE'  3,171 

171  FORMAT ( 37H  NO  FALSE  INTERSTITIALS  WERE  FORMED 
GO  TO  174  . 

172  WRITE  OUTPUT  TAPE  3 , 1 7 3 , M I  ST AKE ( 4 ) 

173  FORMAT (13,  34H  FALSE  INTERSTITIALS  WERE  FORMED  /) 

174  WRITE  OUTPUT  TAPE  3 , 1 75 , MI  ST AK E ( 5 ) 

175  FORMAT ( 10H  WE  FORMED  ,IL,22H  TRIPLE  INTERSTITIALS  ) 

WRITE  OUTPUT  TAPE  3 , 1 08 , I  X , 1 Y , I Z 

v 1080F0RMAT(35H  THE  DIMENSIONS  OF  THE  LATTICE  ARE 


OCCURED  ALONG  AN  ECGE  OR 


ALONG  AN  ECGE  CR  CO 


/  ) 


16H  DIFFERENT  TIMES 


/) 


106 

X 


1 12, 7H  UNITS 
212, 7H  UNITS 
WRITE  OUTPUT 
106  FORM AT ( 87H 
1  Ml 

DO  1025  L  = 

1025  WRITE  OUTPUT 
103  FORMAT! IX, I  4, 7F10. 4, 01  9) 

9998  ENA610),AJP(9990),SLS1(  1 

9990  WRITE  OUTPUT  TAPE  3,9991 

9991  FORMAT!///  32H  THE  RUN  WENT 

9999  JT=- 1 

WRITE  TAPE  6 , JT 

REWIND  6 

END 


/ 1 3H  Y  DIRECTION  12, 7H  UNITS 
) 

TAPE  3, 

L 

E  LB ( L ) /  ) 

1  ,  LMAX 

TAPE  3 , 1 0 3 »  L 


// 1 3H  X  DIRECTION 
/  1  3H  Z  DIRECTION 


VX 


V  Y 


( 6 (  I  , L  )  1=  1 ,7) ,OCTAL(L ) 


16) 


THE  FULL  0 1  STANCE 
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GOX  31 


SUBROUTINE  SITE(V) 

DIMENSION  R ( 7  « 1 5CC ) f  LB ( 1 5CC )  .OCTAL (1  5CC),FI(3),IB!3),T(20) 
DIMENSION  TMINI2C)  ,  I  FACE  (3 ), MI  STAKE!  1C  )  ,MOCN(  15  )  ,NGON(  101.MH1 
DIMENSION  XC ( 6  ) , XP ( 20 ) , Z Z ( 1 5 ) , CON ! 4 ) , C E V ( 7 , 2 ) , A B ( 3 ) 

DIMENSION  DST  ANC  E ( 20 ) 

COMMON  B.LB.OCTAL.FI, I B , T , TM  I N , I  FACE  ,M I S T AK E , MOON .NOON , W , XC  ,  X P  ,  Z Z 
COMMON  CON.DEV, AB, A.NPLANC ,N LINE, LMAX,MASK1,MASK2,MASK3, MAS K4 
COMMGN  I  TEMP, JTE MP, MT EMP ,N TE MP , IN2,IN3,IL,JA,JB,JC,JK,JL,JM,JP,JR 
COMMON  JT,JW,MM,KMIN,MAX,TEMP,NI,JPMAX,NIF,IX.IY,IZ,NN,MAXLAT,JJ 
COMMON  IP,TE,TEM,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,T1,T2 
COMMON  DSTANCE, FMASS, 13, TSL 
COMMON  MASK5,MASK6,MASK7, MAS K8, MOVE,  NCV* 

EQUIVALENCE  ( XP , TM I N ) , { L e , CC T AL ) 

JT  =  M 

LDA( JT), I N  A ( —  1 ) , ENG ( OR ) ,DVI (NPLANE ) ,  STA(IR  +  3),LLS!48),ENC(0) 

DVI (NLINE) ,STA!  IB  +  2  )  ,CLS  (  1 ) ,STQ(  IB  +  1  ) , ADD ( I e  +  3 ) , ENQ! 1 ) , STL ( JT  ) 
LDA( JT) , AJP( 234C ) , R AO ( IB+1 ) 

2340  RETURN 
END 

BOX  32 

SUBROUTINE  INVAC(MJ) 

DIMENSION  B(7, 15C0) ,LB(1500) .OCTAL (1 500) ,FI (3) , IB! 3) ,T(20> 
DIMENSION  TMIN! 20) ,  I  FACE (3) .MISTAKE!  1C) .MOON! 15) .NOON! 10  )  ,W(4  ) 
DIMENSION  XC ( 6  ) , XP ( 2 C ) , Z Z ( 1 5 ) ,CCN ( 4 ) , DE V ! 7 , 2 ) , AB ( 3 ) 

DIMENSION  DSTANCE ( 20 ) 

COMMON  B, LB, OCTAL, FI , I R , T , TM I N , I  FACE , M I S T AK E , MOON , NOON, W , XC.XP.ZZ 
COMMON  CON, DEV, AB, A, NPLANE .NLINE, LMAX.M A SKI , M ASK2 , M A SK 3 , M A S K 4 
COMMON  ITEMP, J TEMP, MT EMP ,N TEMP, I N2, I N3, IL,JA,JB,JC,JK, JL.JM.JP, JR 
COMMON  JT, JW,MM,KMIN,MAX,TEMP,M ,JPMAX,NIF,IX,IY,IZ,NN,MAXLAT,JJ 
COMMON  IP,TE,TEM,C1 ,C2,C2,C4,C5,C6,C7,C8,C9,C10,C1 1 ,C12,T1,T2 
COMMON  DSTANCE, FMASS, 13, TSL 
COMMON  MASK5.MASK6, MASK7 .MASK8.M0VE, NOW 
EQUIVALENCE  ! XP , TMI N) , ! LB, OCT AL ) 

2360  IN2=0 
MM=  M  J 

LILKMM)  ,LDA1  (LB),SCL  (MASK2)  ,SCL  IMASK6)  ,  SST  (  MASK7  ) ',  ST  A 1  (  LB  ) 

2361  NIF  =  0 

LDA {MIST  AKE  +  2 ) ,STA(MISTAKE+10) 

DO  2370  1=1,3 
I F AC  E ! I ) =0 

I F ( AB { I ) — A-B ( I. KM) )  2375,2375,2365 
r  2375  ENA1 !CB) , ALS! 1 ) , ST A  1 { I F ACE ) , RAO ( N IF )  ,SLJ!2370) 

2365  IF! B! I ,MM) — A)  2380,2370,2370 

2380  ENA1 (OB) ,ALS< 1 ) , INA ( - 1 ) ,STA1 {I  FACE), RAC (N  IF  ) 

2370  CONTINUE 

DO  2387  1=1,3 
IB! I  )=8( I ,MM) /A 
TEMP=I B ! I ) 

I F ( B ( I , MM ) / A—  TEM P—0 . 5 )  2  367,2386,238  6 

2386  RA01 ( IB ) 

2387  CONTINUE 
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BOX  33 


IL=IB(3)*NPLANE+IB(2)*NLINE+IB(1)/2+l 

LDA(  IB+1  ) »  ADD (  IP, +2)  ,ADD(  I  P  +  3)  ,ENG(  1)  ,STL(ITEMP)  ,LDA(  ITE”P) 
AJP1(2410),STQ(JM),SLJ(2415) 

2410  JK=0 

2415  LDA ( I B  +  2 ) » ADD ( IB+3),ENC(1),STL(JL),RSC(JL),LCA(JM) ,ADD( JL) 

+  AJP3 { L  +  3 ) * LDA( IL ) *  STA(NCW) , LDA (MOVE ) , A JP 1  ( 274 1  ) , SLJ ( 24  3  2 ) , RAC ( I L ) 
2432  J  A  =  C 

DD  2115  1*1 » 3 

LDAl (  I  FACE) ,AJP( 21 15) , AJP3I2 11 51  ,RA0 ( JA) , INA(-1 ) ,STA( JW) 

AJP1  ( 2120) ,LDA1 ( I FACE) , S TA ( JB ) , ENA ( 1  )  ,  S TA ( JD ) 

2115  CONTINUE 

LDA( JA) , AJP( 2125) ,SLJ(2140) 

2120  RAOINISTAKE+2) »SLJ(214C) 

2125  ENA ( 0 ) , STA ( JB ),STA(JD),LCA(JM),STA(JC),SLJ(214l) 

2140  JC=  1 

214  1  ENA ( 1 ) ,STA( JK ) ,LDA ( JM) ,MLI ( 6B)  ,MUI ( JD ) , I NA(  1 ) , ADD( JB) ,  ADC ( JC ) 

STA(JW),INA(-2),AJP(2447),AJP2(2448), ENA (6) , STA ( JPMA X ) , SL J ( 2 700 ) 
2447  ENA{ 1 5B) ,  STA (JPMAX) , SLJ ( 2700 ) 

24  48  LDA  (  JW)  ,  INA(-1  IB  ),AJP3  (2449)  ,  ENA  (  1  IB  ),  STA(  JPMAX)  ,SLJ  (2700 
2449  ENA(5),STA( JPMAX) 

2700  CO  2698  1=1,10 

2698  NOON ( I )=0 

DO  2699  1=1,15 

2699  M  0  0  N  ( I )=0 

BOX  34 

ENI  1  ( 1  ) 

DO  2475  JP=1  , JPMAX 

GO  TO  (2701,2702,2703,2704,2705,2706,2707,2708,2709,2710,2711, 
12712,2713,2714).  JW 

2701  GO  TO  (2455,2451,2452,2453,2454,2465)  ,JP 

2702  GO  TO  (2455,2451,2470*2456,2452,2453,2457,2458,2459,2454,2460, 
12461 ,2462), JP 

C2500  THIS  SECTION  IS  USED  WHEN  THE  PARTICLE  SLCWS  DOWN  WITHIN  CNE  UNIT 
OF  AN  EDGE.  WE  CAN  NOT  CONSIDER  ALL  THE  NEAREST  NEIGHBORS  IN  A 
CASE  LIKE  THIS,  ONLY  CERTAIN  SELECTED  ONES,  DEPENDING  ON  WHICH 
FACE  IT  IS  THAT  WE  ARE  NEAR. 

2703  GO  TO  (2455,2451  ,2452,2453,2454) ,JP 

2704  GO  TO  (2451 ,2452,2453,2454,2465) ,JP 

2705  GO  TO  (2455,245 1 ,2453,2454,2465) ,JP 

2706  GO  TO  (2455,2452  ,2453,2454 ,2465) ,JP 

2707  GO  TO  (2455,2451 ,2452,2453,2465) ,JP 

2708  GO  TO  (2455, 2451, 2452, 2454, 2465), JP 

2709  GO  TO  (2455,2470,2456,2457,2458,2459 ,2460,2461 ,2462) ,JP 

2710  GO  TO  (2455,2451  ,2452, 2453, 2454, 2458, 2459, 2461, 2462), JP 

2711  GO  TO  (2455,2451 ,2470, 2453, 2457, 2453, 2454, 2460, 2461), JP 

2712  GO  TO  (2455,2452,2453,2454,2456,2457 ,2459,2460,2462) ,JP 

2713  GO  TO  (2455,2451 ,2452,2453,2456,2457 ,2458,2459,2470) ,JP 

2714  GO  TO  (2455,2451 ,2470,2456,2452,2454 ,2460, 2461, 2462), JP 
2451  LDA( JL) ,ADD INLINE ) , ADD ( I L ) ,SLJ(2435) 

24  52  LDAl IL)  , SUB  INLINE) ,ADD( JL) ,SLJ(2435) 

245  3  LDAl  IL) , ADD( NPLANE) , ADD ( JL ) , SL J ( 24 35  ) 

2454  LDA (  I  L  ), SUB ( NPLANE ), ADC ( JL ), SLJ ( 2435  ) 

2455  LDAl I L ) , SLJ ( 2435 ) 

24  56  LDA ( IL) ,SUB( NLINE) . INA( 1 B) , ADD( JL)  ,SLJ(2435  ) 

2457  LDAl IL) , ADD (NPLANE ) ,INA( IB) , ADD( JL) , SLJ (2435 ) 

2458  LDAl  IL), ADD (NPLANE) ,ADD(NLINE)  ,SLJ(2435) 

24  59  LDAl  I  L  )  ,ADD ( NPLANE )  , SUB ( NL I NE  )  , S L J ( 2 435 ) 

2460  LDAl  IL) , SUB (NPLANE) ,INA(  1 B ), ADC ( JL ), SLJ ( 2435  ) 

2461  LDA (  IL)  ,SUB( NPLANE) , ADC( NL I NE ) , SLJ (2 4  35  ) 

2462  LDAl IL) , SUB ( NPL ANE ) , SUB ( NL IN E ) , S L J ( 2 435 ) 

2465  LDAl  ID,  INA(-l)  ,SLJ(2435) 

2470  LDAl  IL)  ,  ADD ( NL I N E )  , ADO ( J L ) , I NA ( 1 ) 

2435  STA( I  TEMP) ,AJP (2475  ), AJP3 (24  75) ,SUB( MAXLAT) , +AJPI L+l ) ,AJP2(2475) 
LDAl  I  TEMP ) ,STA1  ( MOON) , INI  1 ( 1 ) 

2475  CONTINUE 
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BOX  35 


INI  1  ( - 1  )  ,  SILK  JPMAX)  , LD A ( MC V  E )  ,AJP1(  274  1  ) 

2480  EN A (  1  )  ,  S  T  A  (  J  P  ) 

00  2485  K=l, JPMAX 

LDA3 ( MOON) , STA ( I  TEMP ) , L I L 1 ( I  TEMP ) , LO C 1  ( L B ) , LDL ( 4 B ) , A J P 1 ( 24 8 5 ) 

ENQ 1 (0) * LIL1 ( JP)  ,  STQ 1  (  NOCK) ,RAQ(  JP) 

2485  CONTINUE 

LDA( JP) , INA(-1)»STA(JK) 

BOX  36 

LDA ( JK ) ,  INA(-1 ) , AJP( 2490) ,  AJP2 (2510)  , LC A ( I N2  ) , A J P 1  ( 2720 ), SLJ ( 2600 ) 
2490  LDA (NOON+1 ) , STA( IL) 

2492  CALL  SITE  (  ID 

2493  DO  2495  1=1,3 
T  EM  P  =  1 8 (  I  ) 

B (  I  ,MM ) =A*TEMP 
2495  CONTINUE 

ENA(3)  ,STA(NI ) ,L I  LI ( IL) , LIL3(MM) ,LOA ( IL) ,ALS(27) ,STA( I  TEMP) 

LDA3 ( LB ) ,LDQ( MA SK2 )  , S SU (  IT  EM P ) , ST A3( LB ), L0A1 (LB),SST(4B),STA1(LB) 
LDA { IN2), AJP1 (2800 ) 

RETURN 

BOX  37 

2510  DO  2512  K=1 , JK 

LDA3 ( NOON ) ,ST A ( I  TEMP ) , LD A ( C . 0 ) , STA3{  CST ANCE ) 

CALL  SITE ( I T EMP) 

DO  2512  1=1,3 
TE MP= I B (  I  ) 

TEMP=TEMP*A 

2512  DSTANCE ( K ) =DSTANCE ( K )  +  ( TEN P-B ( I , MM) )  * { TEMP-B ( I , MM )  ) 

EN I  2 ( 1  ) 

DO  2514  K  =  1 , J K 

LDA3 (DSTANCE ) , F S B2 ( DST ANCE ) , A J P ( 2 5 1 4 )  ,  A J P2  (  2514),SIL3(ITEMP) 

LIL2( ITEMP) 

2514  CONTINUE 

EN  I  5  (  1  )  ,LDA2(N00N) ,STA5( MOCN ) , INI  5(1  ) 

DO  2516  K= 1 , JK 

LDA2( DSTANCE) ,  F SB3 ( DST ANCE ) ,  +  A JP ( L  +  l  ) , A J P2 ( 2 5 1 6 ) , +S I L2 ( I T EMP ) 

'  LAC(  ITEMP) , I N A3 ( 0 ) , AJ  P ( 2  5 1 6 ) ,LDA3( NOCN) , STAS ( MOON ) , I N I  5  (  1 ) 

2516  CONTINUE 

I N I  5 ( - 1 ),ENA5(0) , IN  A ( -  1 ) ,AJP1 (L+3) ,+ LDA2 ( NOON) , S T A ( I L ) , +SL J ( 2492 ) 
+ENA5( 0) , STA (MMAX) 

DO  2518  K=1 ,  MMAX 

LDA ( MM ) , SUB3 ( MOON ) ,AJP(25?0) 

2518  CONTINUE 

2520  ENA 3 ( 0 ) , AJ  P { L  +  3 ) ,LDA3 ( MOCN ), STA ( IL),+SLJ(24  92)  , +ENQ6 ( 0 ) , LCL ( M ASK 5 ) 
+AJP1 (L+2) ,LDA(M0CN+2) , S TA ( I L ) , SL J ( 2 492 ) ,+LDA(MCCN+l ) , ST A ( I L ) 

SLJ ( 2492 ) 
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BCX  38 


2600 

2605 

2620 


1700 

2621 


2625 


2630 

2635 

2640 

2645 


2650 

2655 

2660 


2665 

2675 

2680 

2686 


T  1  =  0 . 0 

DO  2605  1=4,6 
T 1 =T 1 +8 ( I,MM)«B( I ,MM) 

ENA(  1  ) ,STA( JA) ,LDA(T1 ) ,AJP1 ( 2620) ,ENA(0),STA(JA),SLJ(262C) 
T2=0 . 0 


DO  2645  K= 1 »  J  PMA  X 

LOA 3 ( MOON ) tSTA( I  TEMP)  , ALS( 27 ) , STA( JT  EMP ) 

DC  1700  1=1 , LM  AX 

LD01 (LB) ,LDL(MASK2) ,SU8( JTEMP) ,AJP1 (  17C0) ,LDL( MASK) ) , AJP (2621  ) 
LDA ( TSL 1 , FMU (  1 . OE  +  3 ) , STA3 ( TM I N ) , SL J( 264  5  ) 

CONTINUE 

CALL  SITE(ITEMP) 

DO  2625  1=1 ,3 
T  EMP= 1 8 (  I  ) 

FI ( I )  =  A*T  EMP 
LDA( JA) , AJP 1 (2635) 

DO  2630  1=1,3 

TMI N(K)  =  TMIN( K)  +  (B( I ,MM)-F I (  I )  )* ( B ( I  , MM) -FI (  I ) ) 

GO  TC  2645 
DO  2640  1=1,3 

T2  =  T2+(8( I.MMJ-F I( I ) )*B( 1  +  3,  MM) 

TMIN(K)— T2/T1 
CONTINUE 

ENA ( 1 ) ,STA( KM  IN)  , LDA ( JPMAX ) , I NA ( -1 ) , ST A ( JK ) 

DO  2660  K= 1 , JK 

IF(TMIN(K)-TMIN(K+1 ) ) 266C ,2660,2650 
I F ( TMIN(K+1 J-TMIN(KMIN) )  2655,2660,2660 


KM  I N  =  K+ 1 
CONTINUE 

LI  LI (KMIN) , LDA 1 (MOON) , ST A ( I P ) , LD A ( TS L ) , F MU ( 1 00 . 0 ) , S T A ( TE M ) 

LDA 1 (TMIN) ,+FSB( TEM) ,AJP3(L+2) ,ENI1 ( 0) , SL J { 2685 ) 

CALL  SITE  (IP) 

DO  2665  1=1 ,3 
TEMP= I  8 (  I  ) 

FI (  I ) = A*  TEMP 

T ( I ) = ABSF ( 8 ( I , M M  ) - F  I  (  I  )  1 

LCA(  JB)  ,  INA(  1 ) . ARS( 1 ) ,STA( ITEMP) ,LIL 1 ( ITEMP) ,LDA(C.O)  ,STA1  (  T  ) 
ENA ( 1 ) .STA(MAX) 

DO  2675  1=1,2 

LOA  1  (  T)  ,FSB1  (  T+  1  )  ,  AJP  (2  675  )  ,  AJP2(2-67  5)  ,S  ILl  (  MAX  )  ,RAO(  MAX  )  . 
CONTINUE 

DO  2680  I  =  1  , LMAX 

LDQ1  (LB) , LDL ( MASK2 ) , ARS ( 27 ) , SUB ( I P ) ,  AJP (2685) 

CONTINUE 

STA(MM) , ENA ( 1 ) , S T A ( I N2 ) , S L J  (  236 1  ) 
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2800  LDA(  IP)  *STA  (  IL)  «  LDA(  MJ  )  ,  S  T  A  (  MM  )  ,  E  N  A(  C  )  ,  S  TA  {  I N2  )  ,  S  L  J  (  249  2  ) 

THIS  IS  THE  SECTION  OF  INVAC  THAT  MAKES  THE  INTERSTITIALS.  THE  REST  OF 
THE  SUBROUTINE  IS  CONCERNEC  ONLY  WITH  VACANCIES 
2720  TEMP=B(MAX,MM)-B  (MAX,HJ) 

CALL  SI TE  ( IP ) 

DO  2725  1=1,3 
T2= I B (  I  ) 

B {  I  , MM ) =  A*T2 
2725  B(  I , M J  )  =  B  (  I, MM) 

I  F ( TEMP )  2735,2730,2730 
2730  B(MAX,MM)=B(MAX,MM)  +0.5*A 
B(MAX,MJ )=B (MAX, MJ)— 0.5*A 
GO  TO  2740 

2735  B(MAX,MM)=B(MAX,MM)-0.5*A 
B(MAX,MJ )=B(MAX,MJ)+0.5*A 

2740  CONTINUE 

LDA { IP) ,ALS (27) , STA { I  TEMP ) , LDA { MM ) , A L S < 1 5 ) ,  A CD (  I  TEMP ) , STA(  JA) 

LD0(MASK3) , L I L  3 ( M  J ) , LD A3 ( LB ) , SSU ( J A)  , ST A3 ( LB ) , ENA3 ( 0 ) , ALS ( 1 5 ) 

ADD ( I  TEMP) ,STA( JA) , LI  LI ( MM) , LDA 1 (LB),SSU(JA),STA1(LB),LIL1(IP) 

LDA1 (LB) ,SST(4E),STA1 (LB) ,LUA(MM),SUB(MJ ),AJP1 (2741 ) 

RAO ( MI STAKE+4 ) 

2741  RETURN 
END 


BOX  40 

SUBROUTINE  I NTER ( MO , I  I  ,  S  ) 

DIMENSION  B(7, 1 500) ,LB{ 1500) ,0CTAL(1 500) ,FI (3) , I B< 3) ,T(20) 
DIMENSION  TM IN ( 20 ) , IFACE (3) , M I  STAKE! 10) ,MOON (15), NOON( 10 ) , K(4 ) 
DIMENSION  XC  (  6  )  ,  XP  (  20.) ,  ZZ  (  1  5  )  ,  CCN  ( 4  )  ,  DE  V  (  7 , 2  )  ,  AB  (  3 ) 

DIMENSION  DST  ANC  E ( 2  0 ) 

COMMON  B, LB, OCTAL, FI ,IR, T , TM I N, I  FACE  ,M I  STAKE , MOON, NOON ,W , XC , XP , Z Z 
COMMON  CON, DEV, AB,A .NPLANE ,NLINE ,LMAX, MASKl , MASK2, MASK3 , MASK4 
COMMON  I  TEMP  ,  JTE  MP  ,  MT  E  MP  ,NTE  MP  ,  I  N2 , 1  N3  ,  I  L  ,  J  A  ,  JB  ,  JC  ,  JK  ,  JL  ,  JM  ,  JP  ,  J  R 
COMMON  JT , JW, MM, KM I N, MAX  ,  TEMP , NI ,JPMAX,NIF,IX.IY,IZ,NN.MAXLAT,JJ 
COMMON  IP,TE,TEM,Cl,C2,C2,C4,C5,C6,C7,C8,C9,Cl0,Cn,Cl2,Tl,T2 
COMMON  DSTANCE,FMASS, 13, TSL 
COMMON  MASK5,MASK6,MASK7,MASK0,M0VE,  NOW 
EQUIVALENCE  ( XP , TM I N ) , ( L E , CC TA L ) 

MM  =  MC 

DO  1C  1=1,6 

XP (  I  )=B( I , MM ) -B ( I, I  I ) 

FMU ( 0 . 5 )  STA ( ZZ+2 )  . 

10  XC  (  I ) =ZZ ( 2 ) 

ZZ ( 9)={ XP{4) *XP (4 )+XP( 5 )*XP(5)+XP(6) «XP(6))*1.0E8 
FMU ( C4 )  STA(ZZ+10). 

ZZ ( 1 )=SQRTF( ABSF(ZZ(9) ) ) 

ZZ  ( 10)  =  1 .O/ZZl 10) 

ZZ (  1  1  )=S»S 
ZZ( 7  )  =C1 
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15  ZZ(2)=C2/EXPF(ZZ(7)*C3) 

F  SB { C9 )  FMUfZZ+lO). 

STAIZZ+8) 

ZZ(2)=C3*ZZ( 2)*ZZ(7)*ZZ(  1 C  ) »  ZZ  1 7 ) 

Z  Z ( 5)=ZZ(7)*{ 1 . C-ZZ { 8 ) ) 

ZZ( 2  )  =ZZ ( 7 )  —  ( ZZ(7)»ZZ(5)-ZZ(  11 ) )/(2.0*ZZ(5)+ZZ{2) ) 
LDA  (  ZZ  +  2  )’ 

FSB ( ZZ  +  7 )  A JP2 ( 19) 

LDA ( ZZ+2 ) 

FMUM. OOOOOl)  F  SB  (  Z  Z  +  7  )  . 

A JP2 ( 20 )  LDA ( ZZ  +  2  )  . 

STA ( ZZ+7 )  SL J ( 15) 

19  Z Z  C  2 ) =ZZ ( 7 ) 

BOX  42 

20  LOA(O.O)  STAIZZ+5)  . 

ST  A ( ZZ  +  6 1  LDA ( ZZ+2 )  . 

FMU ( ZZ  +  2 )  S  T  A { ZZ  +  7)  . 

ZZ{ 8)=C2/EXPF{ZZ (2)*C3 )-C9 
LACIZZ+2)  F  M  U ( C  8 ) 

F AD (  1.0)  STA ( ZZ  +  4  )  . 

DO  25  1=1,4 

ZZ  t  9 ) =ZZ ( 4 ) *CON  ( I ) 

ZZ ( 3  )  =  ZZ ( 2 ) »ZZ( 2 ) 

ZZ ( 15)=ZZ(3)*ZZ(10)/ZZ(9) 

ZZ(  1 4 ) =Z Z ( 1 1 ) * ( 2 . 0- ZZ ( 9 )  ) 

ZZ  C  9 )  =  1 .O-ZZ (9) 

ZZ( 1 2 ) = 1 . 0/ ZZ ( 9 ) 

ZZ(9)=C2/EXPF(ZZ (2)*ZZ(12)*C3)-C9 

Z  Z ( 9 ) =W ( I ) /SQRTF { ABSF (ZZ(14)+ZZ(15)*(ZZ(8)-ZZ(9)))) 
FAD ( ZZ  +  5 )  ST  A ( ZZ+5 )  . 

LDAIZZ+9)  FKU(ZZ+12>. 

FMU { Z Z  + 1 2 )  FAD ( ZZ  +  6  )  . 

STA ( ZZ+6 ) 

25  CONTINUE 

ZZ(7)=2.0*SQRTF( ABSF(ZZ(4  )  )  ) 

ZZ(5)=S*ZZ(5)*ZZ(7) 

ZZ( 6)=ZZ(3)»ZZ(6)*ZZ(7)/ZZ(1 ) 

BOX  43 

ZZ(7)=SINF(ZZ(5) )*C6 
ZZ ( 8 ) =COSF ( ZZ { 5 ) ) 

ZZ(3)=ZZ(5)+ZZ(13) 

ZZ( 4)=SINF(ZZ(3) )*ZZ( 1 J+C.5E-4 
ZZ(  1 0 ) =C  5*ZZ ( 1 )*C0SF(ZZ(3) ) 


BOX  44 


F I (  1)=XP(2)*XP(6)-XP(3)*XP(5) 
FI(2)=XP(3)*XP(4)-XP( 1 ) *  X  P { 6 ) 
FI(3)=XP(l)*XP(5)-XP{2)sXP(4) 
ZZ  (  1  )=FI(2)*XP(3)-FI(3)*XP(2) 
ZZ(2  )  =  F  I ( 3 ) *X  P ( 1 > — F I (  1  )  *  XP ( 3 ) 
ZZ( 3)=FI(1 ) *XP ( 2 )-F I ( 2 ) *  XP ( 1  ) 
FMU ( ZZ  +  3 )  STAUEMP) 

LOA ( ZZ+2 )  FMU(ZZ+2) 

FAD(TEMP)  STAUEMP) 

LDA ( ZZ+ 1 )  FMU ( Z  Z  + 1 ) 

FAD(TEMP)  STA(TEKP) 

TEMP=SQRTF( ABSF( TEMP) ) 

L  DA ( 1.0)  FDV(TEMP) 

STA ( TEMP )  LDA ( 2  Z+  1 ) 

FMU ( TEMP )  STA ( ZZ+ 1 ) 

LDA ( ZZ+2 )  FMU ( TEMP ) 

STA ( ZZ  +  2 )  LDA ( Z  Z  +  3 ) 

FMU ( TEMP )  STA ( ZZ+3 ) 


BOX  45 


DO  30  1=  1,3 
13=1+3 

ZZ ( 5  )  =XC ( I)*ZZ(8)+ZZ( I ) *  ZZ ( 7 ) 

ZZ(9)=B( I , I  I )+XC( I >  +  { B { 13,1 I)+XC( 13)  )«ZZ(6) *1.0E14 
DEV (  I, 1 )=ZZ(9)+ZZ(5) 

CEV(  I,2)  =  ZZ(9)-ZZ(5) 

ZZ( 5  )  =  XC( I)  +  ZZ(10)+ZZ(I)«ZZ(4) 

ZZ ( 9  )  =B ( 13, I  I )  +  XC( 13) 

DEV (  13, 1 )=ZZ(9)+ZZ(5) 

30  DEV (  13,2 )  =  ZZ ( 9 )-ZZ ( 5 ) 

LDA(C.O),STA(Tl ) ,STA(T2) 

DO  29  1=4,6 

T 1  =  T  1  +D6V ( I , 1 )*DEV( 1,1) 

29  T2=T2  +  DEV( I ,2 )*DEV{  I,  2) 

DEV (7,1 )=FMASS«T1*0.5180C761500 
DEV(7,2)=FMASS*T2*0.5180C761500 
END 
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BOX  1 


PROGRAM  SLAVE 

DIMENSION  A8(3) ,e<7,2370),FI(3), I B { 3 ) ,L B ( 237C ) , XVZ ( 3 ) , EMAX ( 30 ) 

Cl  MENS  ION  YENTR V ( 1C ) , ZEN7RY (1 0 ) , CC TA L.(2  370 ) .KET ( 26 ) 

DIMENSION  IHIST0(105),NUM8ERS(56).MISTAKE(10)  ,  . 

DIMENSION  NSIDE(6),ESIDE(6) ,N I  TEMP (10) , IC( 1 00 ) . TSL I ( 30 ) , NTSC ( 3C ) 
COMMON  IX, I Y,  IZ,FI,A,NLINE,KPLANE,LMAX,IB,LB,XYZ,MAXLAT, B, OCTAL 
COMMON  IH ISTO, NUMBERS, AB, I L, I  TEMP, ITEM, IT ,TEMP,T£M,KET 
COMMON  NSU8ST ,EBUL»  ETHER M, Ns  IDE ,NBUL , NS  A ME, INTERST.INPAI RS ,NFACE 
COMMON  ESIDE.NV AC, ECAL,OENERGY,N ITEM P, DEPTH, PERC 1ST, IC,EPB 
COMMON  YENTRY , ZENTRY 
EQUIVALENCE (LB, OCTAL) 

REWIND  8 
PRINT  109 

REAO  TAPE  8,  IX,  I Y,  I Z  ,  AB  ,NL  INE  .NPLANE  , MAXLAT ,LMAX,A ,  SPHI ,  FMASS', 

1 ETHRESH, THERMAL, EPB, LULL, ALPHA, BETA,  NTS, MNTS . TT IME, TSLI , NTSC, 
20ENERGY,TENERGY, MISTAKE,  YENTRY,  ZENTRY, NA, EMAX, EFOUNO 

REAO  TAPE  8  , ( LB ( L), ( 8 ( I ,L)  1*1,7)  L*1,LMAX) 

WRITE  OUTPUT  TAPE  3,398 

398  FORMAT { 30H  SET  SLJ  SWITCH  2  UP  IF  NEEDED) 

SLS  C  399  )  . 

399  CONTINUE 
PRINT  101 
DO  400  1=1,6 
READ  INPUT  TAPE  2,100 


100  FORMAT ( 80H 


400 

405 


200 

210 

220 

101 


1 


) 


PRINT  100 

READ  INPUT  TAPE  2, 405, )( KET  ( I )  1*1,25) 

FCRMATMOX,  1314) 

KET { 26 )= 1 000 
SLJ2( 200 )  SLJ (220) 

PRINT  210 

FORMAT ( 40H  SELECTIVE  JUMP  SWITCH  NUM8ER  TWO  WAS  UP) 
PRINT  101 
FORMATt ////) 

IHISTOt 1 )=5H( 13, F 
I H I ST0(2)=4H7. 1 , 

IHI STO( 3 ) =5H3H  P. 

IHIST0(4)=6H4HC. 

IHI STO { 105)= IH ) 

I TEMP=32768 

I2E  1 5  =  I  TEMP 

I2€24=ITEMP«512 

F2E24=I2E24 

I2E30=ITEMP*ITEMP 

MASK25=77770B 

MASK69 =7 777000008 

MASK  13=77770000000008 

MASK14=70000000000000B 

MASKSA=1 OOOOOCCCCOCOOOB 

MASK$U=200000000000000B 

MASK SET =73000000000000078 
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BOX  2 

+  LDA (NTS)  AJP1 ( L+2)  LDA(MNTS)  STA(NTS) 
NSUBST*0 
EBUL=0 • 0 
ETHERM=E FOUND 
DC  1 1105  1*1,6 
ESIDEl  I ) =0.0 
1405  NS  IDE ( I )*0 
NBUL=0 
NSAME^O 
INTERST*0 
NAA=NA— 1 


DO  1407  1-1, NAA 
1407  TSLI(I)*TSLI(I)*1.0E— 14 
TTfME*TTIME*1.0E-14 
DO  1410  L* 1 » LHAX 
+LDA4 ( LB )  SCL(MASKSET)  SlJ2(l+2> 
1410  CONTINUE 


SCL (MASK69 )  SCL (MASKl 3) 


♦STA4ILB) 
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V 


1411 

1412 

1415 

1416 
141? 


1425 


1426 

1427 


1428 

1444 

1445 


1449 


1450 

1451 

1452 


ECX  3 


!:2i 


ICI 

ICI 


( NFA 
(NFA 


W 


6  (  7 »  L  ) 


DO  1449  L=1,LMAX 
LDQ41LB)  LDL ( MASK  1 4 ) 

AJP (1415)  ALS ( 9  ) 

STA(NFACE) 

NSIOE(NFAC( 

ESI DE (NFACc 
GO  TO  1449 
SLJ2( 1416) 

CALL  ROUNO(L) 

+  SL  J ( L+3 ) 

+LDQ4 ( LB )  LDL (MASK13 )  ARS(27)  STA(IL) 

LILKIL)  LDA(IL)  ALS(27)  STA(IL27) 

SLJ  2 (  1417)  LDA4 ( LB )  SST(IL27)  STA4ILB) 
TEMP-BI7,L)-THERMAL 
A JP3 ( 1 425 )  RAC(NBUL) 

EBUL=EBUL+B(7.L) 

LCA4 ( LB )  SST(2B)  STA4ILB) 
TEM=B{7»L)-ETHRESH 
A JP  2 (  1449) 

CALL  ROUND ( L ) 

LIL 1 ( IL) 

LDAl(LB)  SST  ( 4B  )  STAKLB) 

GO  TO  1449 
ETHERM^ETHERM+BITtL) 

SL J2 (  1445  )  LOCI  (LB ) 

LDL ( 4B  )  AJP ( 1 427 ) 

DO  1426  J=1 » L 

LDQ2 ( LB  )  LDL( MASK13) 

SUB (  I L27 )  A J P 1 (  1426) 

LDQ2ILB)  LDL ( 28 ) 

AJP1(  1426) 

ITEMP=L 

AL  S (  15)  ADD21LB) 

I TEMP=J 

ALS (15)  ADD41LB) 

SLJ (1444) 

CONTINUE 
LDAKLB) 

LAC( IL) 

ADD4 ( LB ) 

LDA ( MASKSA) 

LDAKLB) 


SCL(MASKSA)  SCL(MASKSU)  STA21LB) 
SCL(MASKSA)  SCL(MASKSU)  ST  A4 (LB ) 


SST  (4B  )  STAKLB) 

INA4 ( C )  AJ  P (  1 428 )  LCA(MASKSU) 

STA4 ( LB )  SLJ ( 1444 ) 

ACD4 ( LB )  STA4 ( LB  ) 

SST  (  4B  )  STAKLB)  SLJ(1449) 


LDAKLB)  SST  (4B )  STAKLB) 

+LDQ4 ( LB )  LDL ( MASKl 3 )  ARS(27)  STA(IL) 

♦  LILKIL)  LDAKLB)  SST(IB)  STAKLB) 

+  LDL  (  MASK69  )  AJPK1449)  LAC(IL)  INA4(0) 

+  AJPKL  +  3)  LCA4  (  LB  )  SST(NASKSA)  ST  A4  (  LB  )  .  SLJ  (  1 449 ) 
+LDA41LB)  SST ( MASKSU)  STA4 (LB ) 

CONTINUE 

NVAC=0 

DO  1456  L= 1  *  LNAX 

LDQ4  (LB )  LDL  ( 4B  )  AJPK145C)  RAO(NVAC) 

LDQ4 ( LB )  LDL (MASKSA )  AJP(1451)  RAO(NSAFE) 

LDQ4(LB)  LDL ( MASKSU )  AJP(1452) 

LDQ4 ( LB )  LDL ( MASK69 )  AJP11456) 

NVAC*NVAC+MAXLAT-LMAX 


RAC (NSUBST ) 
RAO (I NTERST) 
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ecx  4 


TEMP=ESI DE l 1 ) 

ITEMP^NSIDE ( 1 ) 

00  1457  1=2,6 
ITEMP=ITEMP+NSIDE( I ) 

1457  TEMP=TEMP+ES IDE ( I ) 

ECAL=ETHERM+EEUL+TEMP 

LMAXCAL=  INTERST+NSAME+NSLBST+NBUL+ITEMP 

PRINT  1458,IX*AB ( 1 ) ,  IY*AB(2) ,  I Z , AB (3 ) , MAXLAT , A, SPH I , FMASS , ETHRESH, 
1THERMAL 

1458  FORMAT!  29H  THE  LATTICE  DIMENSIONS  ARE  -5X, 1 1HX-DIRECT ICNI6,  6H 

1  X  A  =F8. 3,  1CH  ANGSTR0MS/34X,  1 1 HY-D IRECT ICNI 6, 6H  X  A  =F8.3,  1  OH 

2  ANGSTR0MS/34X,  11 HZ-DI RECT I  ON  16 , 6H  X  A  =F8.3,  10H  ANGSTROMS//  21H 

3  THE  LATTICE  CONTAINS  1 5 , 6H  AT0MS//33H  THE  LATTICE  SPACING  CONSTANT 

4  A  =F6 .3,  1CH  ANGSTRCMS/27H  THE  SPHERE  OF  INFLUENCE  ISF12.3,  1CH 

5  ANGSTROMS/  19H  THE  ATOMIC  PASS  ISF20.3,4H  AMU// 

6  40H  THRESHOLD  ENERGY  FCR  VACANCY  P RCCUCT ICNF8 . 3 , 4H  EV.  / 

7  25H  THERMAL  THRESHOLD  ENERGYF23. 3, 4H  EV.) 

LDA { ALPHA )  FMU(IBO.O)  FCV ( 3 . 14 15926 536 )  ST  A ( YANGLE  ) 

LDA{ BETA )  FMU(180.0)  FCV ( 3. 1 41 5926536 )  ST  A ( ZANGLE ) 

IT=LMAX-MAXLAT 
PRINT  101 

PRINT  1454, I T, LULL, EPB, ALPHA, YANGLE , BET  A, ZANGLE 

1454  FORMAT ( 23H  CONDITIONS  OF  THIS  RUN///I3,  20H  ATOMS  WERE  SHOT  I  Km 

128H  ONE  ATOM  WAS  SHOT  IN  EACHl3,llH  TIME  STEPS/ 

3  38H  ORIGINAL  ENERGY  CF  EACH  INCOMING  ATCMF12.3,4H  EV./// 

4  28H  DIRECTION  OF  INCOMING  ATOMS/ 

526H  ALPHA  =  ARC  TAN (  VY/VX  )  =F8.4,11H  RADIANS  =F8.3,8H  DEGREES/ 
626H  BETA  =  ARCTANI  VZ/VX  )  =F8.4,11H  RADIANS  *F8.3,8H  OEGRE<€S> 

I TEMP=MAXLAT  + 1 

PRINT  1453, ( L ,YENTRY( L-MAXLAT) ,  ZENTR Y { L-MAXLAT )  L= ITEMP, LM AX ) 

1453  FORMAT ( ///37H  POSITIONS  CF  ENTRY  CF  INCOMING  AT0MS//7H  NUMBER,9X, 

1 1HX9X,1HY9X, 1HZ//(I6,8X,5HC.C0C,2F10 .3) ) 

PRINT  1455,NTSfcMNTS,TTIME, t  TSLI ( I ) ,NTSC( I )  1  =  1, NAA) 

1455  FORMAT { 1  Hi , //  41H  THE  NUMBER  OF  TIME  STEPS  OF  THIS  RUN  WASIlO/ 

1  47H  THE  LIMIT  SET  FOR  THE  NUMBER  OF  TIME  STEPS  WASI 4/  31H  THE 

2T0TAL  TIME  OF  THIS  RUN  W AS  1  PE  1 2. 5 ,8H  SECONDS//  41H  TIME  STEP  LENGT 
3H  NUMBER  OF  FIRST  TIME/4X,  3oHI N  SECONDS  STEP  OF  THIS  L 

4ENGTH/1/1PE14.4.I18)) 

LDA(NTS)  SUB ( 30  )  +AJP2(L>2)  LDA ( NT  S)  STA(IT)  SLJI5541  ) 

LDA ( 30 )  STA(IT) 

5541  CONTINUE 
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PRINT  1459,ECAL,TENERGY,CENERGY,ETHERM, ( I,EMAX<  I)  1=1,  IT) 

1459  FORMAT! 1H1 ,/28H  THE  ENERGY  ACCOUNTED  FCR  ISF14.3.4H  EV./ 

1  28H  SUM  OF  ENERGIES  OF  ATOMS  F14.3.4H  EV./ 

2  28H  ENERGY  SHOT  INTO  LATTICE  F14.3.4H  EV./ 

3  45H  ENERGY  ABSORBED  BY  LATTICE  AS  THERMAL  ENERGYF1 2 .3, 4H  EV.//. 

4  29H  NUMBER  OF  ENERGY  OF  MOST/ 

5  29H  TIME  STEP  ENERGETIC  ATOM// (17 , F 16 .3, 3H  EV)) 

I  NT  ERST= I NTERST/2 

PRING  1460,LMAXCAL,LMAX,NSAME»NSUBST,INTERST,NVAC 

1460  FORMAT! 1H1,/38H  THE  NUMBER  CF  ATOMS  ACCOUNTED  FCR  IS  15/ 

1  38H  THE  TOTAL  NUMBER  OF  ATOMS  IS  15// 

2  38H  THE  NUMBER  IN  ORIGINIAL  POSITION'S  IS  15/ 

3  38H  THE  NUMBER  OF  REPLACEMENTS  IS  15/ 

4  38H  THE  NUMBER  OF  INTERSTITIAL  PAIRS  IS  15/ 

5  38H  THE  NUMBER  OF  VACANT  LATTICE  SITES  ISI5////) 

I F( M I  STAKE (  1 ) )  9955,9955,9956 

9955  PRINT  9957 

9957  FORMAT ! 30H  THE  ENERGY  CHECK  DIC  NOT  FAIL  ) 

GO  TO  9965 

9956  PRINT  9959, M I  STAKE ( 1 ) 

9959  FORMAT! 24H  THE  ENERGY  CHECK  FAILED  I4,16H  DIFFERENT  TIMES) 

9965  IF { MI  STAKE! 2 ) )  9961,9961,9962 

9961  PRINT  9963 

9963  FORMAT { 62H  NO  INTERSTITIALS  OR  VACANCIES  OCCUREO  ALONG  AN  EDGE  OR 
1CCRNER  > 

GO  TO  9966 

9962  PRINT  9960, M I  STAKE t 2 ) 

9960  FORMAT  ( 60H  AN  INTERSTITIAL  CR  VACAICY  CCCURRED  ALONG  AN  EDGE  OR  CO 
•  1RNER  I  4 , 1 6H  DIFFERENT  TIMES) 

9966  CONTINUE 
PRINT  101 
PRINT  1462 

1462  FORMAT!  24  X,  6HNUMBER  ,5X',12HT0TAL  ENERGY/5H  FACE,6X,5HGR0UP,7X, 

18HIN  GROUP, 7X,6HIN  EV.//) 

PRINT  1464,.NBUL,EBUL, (NSIDE! I ) ,ESIDE 1 1 )  1  =  1,6) 

1464  FORMAT!  19H  KNOCK-CNS  I9,F16.  3/19H  1  FRONT  FACEIS',F16 

1.3/19H  2  BACK  FACE  I9,F16.3/19H  3  LEFT  FACE  I9,F16;3/ 

2  19H  4  RIGHT  F  ACEI 9  ,F1  6.  3/19H  5  TOP  FACE  I9,F16.3/ 

3  20H  6  BOTTOM  FACE  18 ,F 16.3) 
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BOX  5 


1465 


•R  OF 
:/) 


PRINT  1 465,NS AME 
FORMAT! lHlt1XtS9H  THE  NUMB 
1SITJCNS  ISI5//10H  THEY  AR 
ITEM- 7777 777000077777  B 
IT-77700007777777776 
J  1  1 
K  =  1 
L=  1 

DO  1466  1=1.10 

1466  NUMBERS ( I ) =0 

DC  1480  M-l.LMAX 
I  Ft  J— 11  )  1474,1467,1467 

1467  PRINT  1 472*1  NUMBERS ( I )  1*1,10) 

DO  1468  1=1,10 

14r68  NUMBERS!  I  )=0 
J=1 
L=L  + 1 

IFTK-501)  1470,1469,1469 

1469  PRINT  109 
K=  1 

L  =  1 

GO  TO  1474 

1470  IF(L-ll)  1474,1471,1471 

1471  PRINT  110 
L=  1 

1472  FORMAT! 10110) 

1474  LDA5  ( LB )  SCL(IT)  ARS!27)  STAUTEMP) 

LAC  1 1  TEMP)  INA510)  AJPK  1476)  LDA5(LB) 
SCL  (  ITEM)  AJPK1476) 

NUMBERS!  J)=.M 
1476  J=J+1 
1480  K=K+ 1 

IN  I  2  (- 1 )  SIL2  (  IT)  LDAUT)  AJP1480) 

480  PRINT  1472, {NUMBERS (I )  1*1, IT) 

BOX  6 


ATOMS  THAT  ARE  IN  THEIR  INITIAL  PO 


PRINT  1560, NVAC 

1560  FORMATdHl, 1X,37HTHE  NUMBER  OF  VACANT  LATTICE  SITES  ISI6///10H  THE 
1Y  ARE  //) 

IT-1 

DO  1565  L=1 , MAXLAT 
L0Q4  (  LB )  LDL14B)  AJPK  1565) 

NUMBERS! IT )=L 

RAOIIT)  SUB  (51)  AJPK1565) 

PRINT  1472, (NUMBERS! IT)  IT=1,50) 

IT-1 

CONTINUE 

PRINT  1472,(NUMEERS(I )  1*1, IT) 


1562 


1565 

1566 
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BOX  7 


1499 

1500 

1502 


1515 

1517 

1518 


1520 


PRINT  109 
PRINT  1499 

FORMAT ( 1 3H  REPLACEMENTS//) 

PRINT  1  500,NSUBST 

FORMAT!  IX, 30HTHE  NUMBER  OF  REPLACEMENTS 
1  INITIAL  FINAL  FENETRAT ION 

DO  1502  1=1,100  • 

IC< I ) =0 

DO  1517  L= 1 »  LMA  X 

LDQ4ILB)  LDL ( MASKSU)  AJP(1517)  LDL IMASKl  3 ) 
CALL  INITIAL(L) 

CALL  OEPTHP(L) 


IS  16/// 
PENETRATION  IN/ 


ARS ( 27 )  STA(IL) 


PRINT  1 51 5,L, IL»  DEPTH.  PERDIST 
FORMAT! 17,1 10, FI 6. 3, FI  8. 3) 


1521 

1525 

1527 


1530 


CONTINUE 
PRINT  109 
PRINT  1499 
PRINT  1520 

FORMAT 1 40H  DISTRIBUTION  CF  X-OIRECTICN  PENETRATICN//2 IX, 

117HR  A  N  G  E18X,6HNUMBER/) 

DO  1521  1=1,30 
J  =  I-7 
TEM- J 

TEM=TEM-0.5 
TEMP  =  TEM+;i.O 

PRINT  1525, J. TEM, TEMP, IC(IJ 

FORMAT! 13, F11.1.9H  X  A  -  16HX  PENETRATION  -F6.1,4H  X  AI11) 
DO  1527  1=1,30 
NUMBERS!  I  )=>IC!I) 

IL=-7 

CALL  HISTOI30) 

PRINT  109 
PRINT  1499 
PRINT  1530 
FORMAT l 35H 
117HR  A 


DO 

IT' 


1531 

1-1 


DISTRIBUTION  CF  RADIAL  PE N6TR ATI 0N//2 1 X, 
N  _  G  E22X.6HNUMBER/) 


1=1,20 


1531  PRINT  1532,1 , IT, I, IC< 1+30) 

1532  FORMAT! 13,  19, 9H  X  A  -  21HRA0I AL  PENETRATION  -14, 4H  X  AI12) 
DO  1533  1=1,20 

1533  NUMBERS!  I  )=IC! 1  +  30) 

CALL  H ISTO { 20 ) 


S5H 

56H 


91 


BOX  8 


PRINT  109 
PRINT  1534 

1534  FORMAT ( 1 9H  INTERSTITIAL  PAIRS//) 

PRINT  1535,INTERST 

1535  FORMAT! 10H  THERE  AREI4.19H  INTERSTITIAL 
1  THE  PAIR  NUMBER  PENETRATION  OF 
2PENETRAT ION  OF  LL 


3  L  LL  OF 

4  IN  X-DIRECT ION 
DO  1538  1*1,100 

1538  ‘ICC  I )  =0 

DO  1550  L=;  1  ,  LMAX 
LDQ4 ( LB  )  LDL (2B  ) 


PA  IRS// 1 1H  NUMBERS  0F/TC4H 
L  PENETRATION  OF  L 

PENETRATION  OF  LL/  1GSH 

SITE.  IN  X-DIRECTICN  IN  RADIAL  OIRECTION 
IN  RADIAL  OIRECTION//) 


AJP 1  (  1550) 


LDQ41LB)  LDL ( MASK69)  AJP (1550)  ARS115) 

ST  A ( LL)  LDL ( MASK  13 )  ARS!27)  STA(IL) 
LILKLL)  LDQltLB)  L0L12B)  AJPK155C) 
CALL  INITIAL ( L) 

CALL  DEPTHP(L) 

TE  =DEPTH 
TEM  =  PERD I  ST 
CALL  INITIAL(LL) 

CALL  DEPTHP ( LL ) 

PRINT  1545,L,  LL,  IL . TE , TEM , OEPTH, PERD I  ST 
1545  FORMAT! 15,16, 19, F16.4,F2C.4,F22.4,F2 1 .4 ) 
1550  CONTINUE 

DO  1555  1=1, ICO 

1555  ICC I)  =  IC(I  )/2 
PRINT  109 
PRINT  1534 
PRINT  1520 

DO  1556  1=1,30 
J  =  I  -  7 
TEM  =  J 

TEM=TEM— 0.5 
TEMP=TEM+:i.O 

1556  PRINT  1525, J, TEM, TEMP, IC (I ) 

CO  1557  1=1,30 

1557  NUMBERS! I )  =  IC( I  ) 

IL=-7 

CALL  HIST0130) 

PRINT  109 
PRINT  1534 
PRINT  1530 
DO  1558  1=1,20 
I T= I— 1 

1558  PRINT  1 532, I,IT,I,IC( I+3C ) 

DO  1559  1=1,20 

1559  NUMBERS! I)=IC! 1+30) 

CALL  HISTOI20) 
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BOX  9 


L7Xr1 HX8X, 1HY8X, 1HZ8X , 2HVX7X, 2HVY7X 
11HX  RADIAL//) 


LDA(NBUL)  A JP( 9951 ) 

1571  PRINT  109 
PRINT  1572 

1572  FORMATl 10H  KNOCK-ONS//) 

PRINT  573 » N8UL 

573  FORMATl 10H  THERE  AREI5.1CH  KNOCK-ONS//) 

PRINT  574 

574  FORMAT { 91 X, 11 HPENETRATI0N/4H 
1 ,2HVZ5X,6HENERGY7X,5HLB( D1CX, 

1576  FORMATl 15, 7F9. 3 , 3X,0 1 5 ,2F8.2 ) 

DO  1577  1=1,100 

1577  ICCI)=0 
IT=1 

00  1578  L=1 , LMAX 
LDQ41LB)  LDL128)  A J P t  1578) 

LDA(IT)  SUB ( 5  1 )  AJP111567) 

PRINT  109 
PRINT  1572 
PRINT  574 
I  T=  1 

1567  CALL  INITIAL (L) 

CALL  DEPTHP(L) 

CALL  ECOUNT(L) 

PRINT  1576, L, (B( I,L)  1=1 ,7 ) .OCTAL (L) .DEPTH, PERDI ST 
RAOl IT) 

1578  CONTINUE 
PRINT  109 
PRINT  1572 
PRINT  1520 

DO  1569  1=1,30 
J  =  I -7 
TEM=  J 

TEM=TEM-0. 5 
TEMP=TEM+1 .0 

1569  PRINT  1  525».J»  T£M,TEMP,  IC  ( I  )- 
DO  1579  1=1,30 

1579  NUMBERS ( I )=IC ( I ) 

I  L*-7 

CALL  HISTOl 30 ) 

PRINT  109 
PRINT  1572 
'  PRINT  1530 

DO  1574  1=1,20 
I T= I— 1 

1573  NUMBERSl I ) = IC ( 1+30) 

1574  PRINT  1532*1, IT, I, ICl I+3C) 

IL“0 

CALL  HI STO (20 ) 

PRINT  109 
PRINT  1572 
PRINT  1580 

1580  FORMAT (23H  DISTRIBUTION  CF  ENERGY// 

115X17HR  A  N  G  El 4X, 6HNUM8ER/ ) 

DO  1680  1=1,25 

1680  PRINT  1 526*1 ,KET( I  )fKET( 1  +  1) »IC( 1+50 ) 

1526  FORMATl 13,17, 21H  EV  -  ENERGY  -16, 3H  EVI1C) 
CO  1681  1=1,25 

1681  NUMBERSl  I  MCII  +  SO) 

CALL  HISTOl 25 ) 
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ecx  10 


,  LMAX 

LDL IMASK14) 
S IL2 ! IT) 

ARS ( 3 ) 


A JP ( 1 595 ) 
SUB ( IT) 
STA ( NTS) 


ALS ( 9 ) 
AJP1( 1595) 


9951  DO  1599  J  =  1 ,6 

IF  IJ-l)  1581,1581,1583 

1581  PRINT  1 582»NS IDE!1) 

1582  FORMAT! 1 H 1 , IX , 1 6HSPUTTERED  ATOMS-I 10 , 22H  ATOMS  WERE  SPUTTERED.) 

GO  TO  1585 

1583  PRINT  1584,J,NSI DEI J) 

1584  FORMAT!  1H1,  IX, 45HTHE  NUMBER  OF  ATOMS  WHICH  EXITEO  THROUGH  SIDE-12,3 
1H  ISIS) 

1585  LDA21NSIDE)  AJP!1599)  AJP311599) 

PRINT  1586 

1586  FORMAT! ///91X,  23HPENETRAT  ION  TIME  STEP/4H  L7X, 1HX8X, 1HY8X , l 

1HZ8X  *  2HVX7X ,2HVY7X,2HVZ5  X,6HENERGY7X,5HLB(L) 10X,22HX  RADIAL 

20F  EXIT//) 

1587  FORMAT! 15 ,7F9 .3 , 3X, 01 5 ,2F8.2 , 19) 

NFACE=0 

DO  1589  1=1,100 
1589  TC( I ) =0 
NTS  =  0 

DO  1595  L=1 
LDQ4 ( LB ) 

STA(NFACE) 

LDL  IMASK25) 

1594  CALL  INITIAL !L) 

CALL  DEPTHP(L) 

CALL  ECOUNT(L) 

PRINT  1587, L, (£( I,L) 

1595  CONTINUE 
PRINT  1695, J 

1695  FORMAT! 1H1,12H  EXITED 
PRINT  1520 
DC  1696  1=1, 3C 
L  =  I -7 
TEM  =  L 

TEM=TEM-0. 5 
TEMP=TEM+ 1.0 

1596  NUMBERS!  I)  =  IC( I ) 

1696  PRINT  1 525, L , TEM  , TEMP , IC  1 1) 

.  I  L=-7 

CALL  HI STO! 30 ) 

PRINT  1695, J 
PRINT  1530 
DO  1697  1=1, 20 
I T= I— 1 

1597  NUMBERS!  I  )  =, I C ( 1  +  30) 

1697  PRINT  1532, I, IT,  I  , I C ( I +  3C  ) 

I  L'  =  0 


1  =  1 ,7), OCTAL! L) ,OEPTH, PERO I  ST, NTS 


FACE  14 // ) 


CALL  H I STO ( 20  ) 

PRINT  1695, J 
PRINT  1580 
DO  1698  1=1 ,25 

1598  NUMBERS!I)=IC(I+50) 

1698  PRINT  1526,I,KET!I),KET!  I  +  1),IC!I  +  50) 
CALL  HIST0!25) 

1599  CONTINUE 
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ecx  li 


107 

108 
109 

m 
112 
1  13 
1  14 

115 

116 

117 

118 

119 

1 

120 


SL J 1(120) 

PRINT  119 
L  =  1 

PRINT  118.L, (8< I,L>  I**  1 , 7 )  .OCTAL  (  L  ) 
FORMAT ( 1H1) 

FORMAT (/) 

DO  117  Li2»LMAX 

IFC  1 0 -» (  {  L —  1  1/101-L+l  )  117,113,113 
IF( 50*( (L-l )/5C)-L+l)  116,114,114 
PRINT  119 
GO  TO  117 
PRINT  110 

PRINT  118,L,(E(I ,L)  1=1 ,7),CCTAL(L) 
FORMAT (  IS.7F9.3,3X,015) 

FORMAT (1H1»///13H  LATTICE  DUMP///4H 
,2HVY7X,2HVZ5X,6HENERGY  7X,5HL8(L)/J 
ENO 


L7X, 1HX8X, 1HY8X, 1HZ8X,2HVX7X 


BOX  12 

SUBROUTINE  ROUNO(L) 

CIMENSION  A B ( 3 ) »B(7 ,2370  ),FI (3) , 18 (3 ) , LB ( 2370 ) , XYZ ( 3 ) , EMAX ( 30 ) 
DIMENSION  YENTRY ( 10), Z ENTRY!  10  )  ,OCTAL(  2370)  ,KET(26) 

DIMENSION  IHIST0C105),NUMBERS(50),MI STAKE (10) 

DIMENSION  NSIDE16) ,ESIDE (6) ,NI TEMPI! 0) ,IC(100),TSLI (30),NTSC( 20) 
CCMMCN  IX, I Y, IZ,FI , A ,N L I NE , NPLANE , LM AX , I  8 ,L 8 ,XYZ , MAX LAT, 8, OCT AL 
COMMON  I HISTO, NUMBERS, A8,I L, I  TEMP, ITEM, IT, TEMP, T EM, KET 
COMMON  NSU8ST , E8UL, ETHERM ,NS I DE , N8UL ,NSAME, INTERST, INPAIRS ,N FACE 
COMMON  ES IDE, NV AC, EC AL,OENERGY,N ITEM P, DEPTH, PERCI ST, IC.EP8 
COMMON  YENTRY ,ZENTRY 
EQUIVALENCE (LB, OCTAL) 

DO  1320  1=1.3 
1320  IBI I )=B( I,L)/A+C.5 

ITEMP=IB( 1 )+I8(2)+I8(3) 

IF.(2*(  ITEMP/2J-ITEMP)  1  330 , 1325,  1  330 
1325  IL=IB(3)*NPLANE+I8(2)*NL INE+I8! 1 )/2+ 1 
RETURN 

1330  NUMBERS ( 1 ) ^  18(3)  *NPL ANE+  18(2)  *NL INE  + ( IB ( 1 ) +  1 ) /2+ 1 

NUMBERS (2 )=  I E ( 3  )  *NPL  ANE ♦  18(2)  *NL INE+ ( I B ( 1 )-l ) /2+ 1 

NUMBERS ( 3 )=  18(3)  *NPL ANE  +  ( I  8 ( 2 )  +  l  ) «NL INE+  18(1)  /2+1 

NUMBERS ( 4 ) =  IB(3)  *NPLANE+ ( I B ( 2 )— 1  )*NLINE+  IB ( 1  )  /2+1 

NUMBERS! 5)=(I8(3)+1 )*NPLANE+  18(2)  »NL INE+  IBM)  /2+1 

NUMBERS(6)=(I8(3)-1 )*NPLANE+  18(2)  *NL INE+  18(1)  /2+1 

TEMP= 1  .0E50 
DO  1340  1=1,6 
IT-NUMBERS! I  ) 

AJP311340)  SUB ( M AXL AT ) 

A JP {  1333)  AJP2( 1340) 

1333  CALL  INITIAL(IT) 

TEM=0. 0 

DO  1335  ITEMP^l , 3 

Fit ITEMP )=XYZ( ITEMP )— BC I TEMP,L ) 

1335  TEM=TEM+FI( ITEMP) *FI ( ITEMP) 

IF! TEM-TEMP )  1337,1340,1340 
1337  IL=IT 

TEMP=TEM 
1340  CONTINUE 
END 
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BCX  13 


SUBROUTINE  INITIAL (L) 

DIMENSION  AB(3) ,  8 ( 7 , 2 370  ) , F I (3), I  8(3) , L8 (237C )  ,XYZ ( 3 ) 
DIMENSION  YENTRYI 10) , ZEN  TRY! 1 0 ) , CCTA L ( 23 70 ) , KET ( 26 ) 
DIMENSION  IHIST0I105) 

-  NSI0EI6) 


EMAX ( 3C ) 


DIMENSION 
COMMON 
COMMON 
COMMON 
COMMON 
COMMON 


05) , NUMBER  SI  50) .MI  STAKE! 10) 
, ESIOE  (6 )  ,  NITEMP  (10)  ,IC(10C 


)  ,TSLI( 30) , NTSCI3C) 


IX,  I  Y,  I  Z,  FI » A  ,NLI  NE»N'PLANE»  LM AX , I  B, LB, XYZ » MAX LA T,  8,  OCTAL 
I H IS TO, NUMBERS, AB ,IL,ITEMP,ITEM,IT, TEMP, T EM ,KET 
NSUBST ,EBUL, ETHER M , NS  IDE , NBUL , NSAME, I  NT ERST, INPAIRS, NF ACE 
ESIDE,NVAC,ECAL,OENER GY, NIT EM P, DEPTH, PER CIST, IC,EPB 
YENTRY , ZENTRY 
EQU  I  VALENCE (LB,OCTAL) 

IFIL-MAXLAT)  1003,1003,1023 
1003  CONTINUE 

IBt3)=(L— 1)/NPLANE 
STQ ( ITEMP  ). 

IB(2)=ITEMP/NLINE 
STQ( ITEMP). 

1 8 { 1 )=ITEMP«2 
ITEMP=IB(2)+IB(3) 

IFC  2* ( I  TEMP /2)- 1  TEMP)  1 0C5 , 1 0 1 0 , 1 0 10 
1005  'I B 1  1  )=IB( 1 )  +  1 
1010  DO  1 C 1 5  1=1,3 
F I  (  I )  =  I B  (  I  ) 

1015  XYZ (  I  )  =A*FI ( l ) 

RETURN 

1023  XYZ ( 1 )=0.0 

XYZ (2)=YENTRY(L-MAXLAT) 

XYZ ( 3  )  =  ZENTRY ( L-MAXLAT ) 

END 


BOX  14 

SUBROUTINE  ECCUNTIL) 

DIMENSION  AB ( 3 ) , B ( 7 .2 370  ) , F I (3), 1 8(3 ) , LB ( 2370 ) , XYZ ( 3  )  , EM AX ( 30 ) 
DIMENSION  YENTRY( 10), Z ENTRY! 1 0 ) , OCTA L ( 2370 ) , KET (26) 

Cl  MENS  ION  IH I STO ( 105) , NUMBER  S ( 50 ) , MI  STAKE ( 10) 

DIMENSION  NS  ICE(6) ,ESICE(6) ,NITEMP(1 C) , I C ( 1 00 ) , TSL I ( 30 ) , NTSC ( 30) 
COMMON  IX,IY,IZ,FI, A, NL I NE , N PL  AN E, LM AX , I B ,LB , XYZ , M AXL AT, B , OCT AL 
COMMON  IHISTO, NUMBERS, AB,IL, I T EMP, IT  EM , I T , TEMP , T EM ,K ET 
COMMON  NSUBST,EBUL,ETHERM,NSIDE,NBUL , NSAME, I NTERST ,INPA IRS  * NFACE 
COMMON  ESIDE, NVAC,ECAL»OENERGY, NITEMP, DEPTH, PEROIST, IC,EPB 
COMMON  YENTRY, ZENTRY 
EQUIVALENCE (LB, OCTAL) 

TEMP=KET (25 ) 

TEMP=B(7»L )— T  EMP 
EN I  1(25)  AJP2(25) 

DO  20  1=1,24 
TEMP=KET ( 1+ 1 ) 

TEMP=TEMP-B(7,L) 

A JP2 ( 25 ) 

20  CONTINUE 
25  RA0KIC  +  50) 

END 
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ecx  15 


1625 

1630 

1631 

1635 

1640 

1645 

1650 


SUBROUTINE  DEPTHP(L) 

DIMENSION  AB(3)»B(7,2370),FI(3),IB(3)»LB(237C),XYZ(3)»EMAX(3C) 
YENTRY ( 10) .Z ENTRY! 1 0 ) , OCTA L ( 23 70 ) 2 KET ( 26 ) 

IH I STO (105), NUMBERS (50), MI  STAKE ( 1  0 ) 

NSIDE(6),ESICE(6),NITEMP(10)» IC( 100 ) ,TSL I (30 ) , NTSC! 30 ) 
IY  ,  IZ.FI ,A,NLINE,NPLANE,LMAX, I B , L B , XYZ , M AX LAT, fi, OCTAL 


DIMENSION 
DIMENSION 
DIMENSION 
COMMON  IX 
COMMON 
COMMON 
COMMON 
COMMON 


IH  I STO,  NUMBERS  »AB  ,  IL,  ITEMP,  IT  EM, IT , TEMP.T EM , KET 
NSUBST  »EBUL»ETHERM,NSIDE,NBUL ,NSAME, INTERST, INPAIRS, NFACE 
ESIDE»NVAC,ECAL,OENERGY,NITEMP,DEPTH,PERDIST, IC,EPB 

_  YENTRY, ZENTRY 

EQU  I  VALENCE ( LB  »OCTAL ) 

DEP  TH=B  ( 1 , L )-XYZ  <  1  ) 

1  =  1 

TEM  P= I 

I  FI  (TEMP-6. S)«A-DEPTH)  1630,  1635,1635 
IF.(  1-29 )  1631,1631,1635 
1  =  1  +  1 

GO  TO  1625 
iIC(  I  ).=  IC  ( I  )  +  1 

PER0IST=SQRTF<(B(2,L)-XYZ(2) )**2+(B( 3,L)-XYZ(3)  )**2) 

1  =  1 

T  EMP= I 

IFfPERDIST-TEMP«A)  1650,  1645,1645 
1  =  1  +  1 

GO  TO  1640 

IC  t  I  +30 )  =  IC ( I *30  >  +  1 

END 


EOX  16 

SUBROUTINE  HI STO (NGROUPS ) 

DIMENSION  AB(3) , 6(7,2370 ),FI(3),I8(3), LB (2370) , XYZ(3) , EM AX ( 30 
DIMENSION  YENTRY (10), ZENTRY ( 1 0 ), OCTAL ( 2370 ) ,KET (26 ) 

DIMENSION  I  HI STO ( 105) ♦ NUMBERS (50 ), MI STAKE(IO) 

DIMENSION  NS  I DE ( 6 ) , ES I DE (6 ) , NITEMP( 1 C ) , I C ( 1  CO ) , TSL I ( 30 ) , NTSC ( 3C ) 
COMMON  IX, I Y, I Z, FI ,A , NL I NE , N PL ANE ,LM A X , I B, LB ,XYZ , MAX L AT, E , OCT AL 
COMMON  IHIST0,NUM8ERS„AB, I L, ITEMP, ITEM, I T ,T EMP , TEM, KET 
COMMON  NSUBST  ,EBUL,ETHERM,NSI DE  »N8UL , NS  A ME, INTERST, INP AI RS , NFACE 
COMMON  ES IDE, NV AC, ECAL,0 ENERGY, NIT EMP, DEPTH, PERCI ST, IC,EP8 
•  COMMON  YENTRY, ZENTRY 
EQUIVALENCE  ( LB',  OCTAL) 

KT0TAL=0 
NUMB  I G=0 

DO  1110  1=1, NGROUPS 
NTOTAL=iNTOTAL+NUMBERS(  I  ) 

IF( NUMBERS ( I )-NUMBIG)  1 1  10 , 1 1 1 0 , 1 1 08 

1108  NUMB  I G=NUMBERS  (  I  ) 

1110  CONTINUE 

PRINT  1109, NTOT  AL 

1109  FORM AT ( / / 1 0X , 1 7H 1 00  PERCENT  ME ANSI  5/ ) 

TCTALN=NTOT  AL 

B I GNUM=NUNB I G 
DO  1113  J=1, NGROUPS 
FNUMBER=N UMBERS  (  J) 

NUMBER= 1 00 . C+FNUMBER/B IGNUM+O .5 
TEMP  =  1  00. 0«-FNUMBER  / TOTAL N 
DO  1111  1=5,104 

1111  ‘I H I S T C  ( I  )  =  3H1  H 
ITEMP=NUMBER+4 

DO  1112  1=5, ITEMP 

1112  IHISTOt I )=3H1HX 
IT= J+I L 

1113  PRINT  IHISTO.IT, TEMP 
END 

END 


97 


PROGRAM  RON 

DIMENSION  AB(3) ,8(3,2370 ),NT (2000 ),P (3) ,X!3,2000) ,S 
DIMENSION  18:3) , FI  13) ,XYZ (3 ), NUMBERS ( 6 ) , YENTRY ( 10), 
+  SLJKL  +  1)  SL J  (  4 ) 

IT**—  1 

WRITE  TAPE  6,  IT,  IT, IT, IT 
4  REWIND  6 
PRINT  3 
3  FORMAT! 1H1 ) 

PRINT  1 

1  FORMAT!////) 

DO  5  1=1,6 
READ  INPUT  TAPE  2,2 
‘  FORMAT ( 80H 


S(  3) 

ZENTRY! 10) 


1 


) 


5  PRINT  2 
PRINT  1 
SLJ2I7) 

RE AO  TAPE  6, IX,  IY, IZ, AB, MAXLAT , A , SPH I , FM ASS , ETHRESH , THER MA L , MNP 
READ  TAPE  6,((B(ItJ)  1=1,3)  J= 1 ,MAXL AT ) , NLI NE, NPLANE 
PRINT  1098,  IX, AB! 1 ) , I Y ,A6 ( 2 ) , I Z , AB ( 3 ) , MAXLAT , A,SPHI , FMASS, ETHRESH, 
1  THERMAL 

1098  FORMAT { 29H  THE  LATTICE  DIMENSIONS  ARE  -5X  ,  1 1 HX-D  IRECT I  ON  I  6 ,  6H 


2 

3 

4 

5 

6 
7 


_  _ 1  UIN  A  U  , 

X  A  =  F  1 1.7,1  OH  ANGSTR0MS/34X,  1  1HY- C I RECT I  ON  1 6 . 6H  X  A  =F11.7,  1  CH 
ANGSTROMS/34 X,  1  1HZ-D IRECT I CNI6 » 6H  X  A  =F11.7,l0H  ANGSTROMS//  21H 
THE  LATTICE  CONTAN I S 1 5, 6H  ATCMS//33H  THE  LATTICE  SPACING  CONSTANT 
A  =F 1 0 • 7, 1CH  ANGSTR0MS/27H  THE  SPHERE  OF  INFLUENCE  ISF16.7,  10H 

ANGSTROMS/  1 9H  THE  ATOMIC  MASS  ISF24.7,4H  AMU// 

4  CH  THRESHOLD  ENERGY  FOR  VACANCY  P RCDLCT I  CNF  1 0 . 3 ,4H  EV./ 

25H  THERMAL  THRESHOLD  ENERGYF25. 3, 4H  EV.) 

PRINT  1 
LMAX=MAXLAT+MNP 

7  READ  INPUT  TAPE  2,8,T 

8  FORMAT(20X,F10.6) 

DX=T»A 

RA=  1 . 0/A 
ENI6( IB) 

10  READ  TAPE  6,IT,P 
LDA(IT)  AJP3! 30) 

CO  15  1=1,3 
T  =  ABSF(P<  I)-B('I,  IL)) 

FSB(DX)  AJP2  (20 
15  CONTINUE 
SLJ (10) 

20  DO  25  1=1*3 
B( I ,  IL)=P( I ) 

25  X ( I ,N)=P( I) 

LDA(IT)  ST A6 ( NT )  +ISK6(3000)  SLJ(IC)  ENI6(3001  ) 

30  'IN  1 6 (- 1 )  ,  S I  L6( I  TEMP )  ENI2(0B) 

DO  45  L=1 ,LMAX 

.......  .  ST  A  (  NO) 


LIU5  ( IT )  SIL5UN)  LIL4UT)  SIL4(IL) 


ENA(-rl  ) 


LIL1 (  IU) 
SIL1 (IN) 


SI L 1 ( I  I  ) 
LCA(NO) 


LDA (  1 1) 
AJP3! 32) 


SUB ( IL  ) 

♦  I SK5 ( 4 )  SL J ( 4  2 ) 


E N 1 5  ( 0 B  )  S  I L4  (  I  L  ) 

CO  4  5  K= 1,1  TEMP 
LDA3 ( NT )  STA(IU) 

AJPK45)  LIUKIU) 

PRINT  41 

41  FORMAT!/) 

LDA ( NO )  A JP2 ( 42 ) 

32  PRINT  35, IL 

35  FORMAT ( // IX  ,20H0RBIT  OF  ATOM  NUMB  ERI 5// 

1 1  OH  NUMBER  0F9X , 1 5HACTUA  L  PC S I T I  ON 16 X , 9HNUMB ER  0F7X , 1 9HNCRMAL I  ZED 
3POSITICN/10H  TIME  STEP6X , 1 HX9X , 1 HY9X , 1  HZ  1 3X , 9HT I  ME  STEP6X, 1 HX9X, 
4lHY9X,lHZ/) 

IN  1 2  (  6) 

N0=  1 

42  DO  40  1=1,3 

40  S ( I )=RA*X( I , K ) 

PRINT  43.IN, ( Xt I ,K)  1=1 , 3 ) , IN , ( S ( I )  1=1,3) 

43  FORMAT! 1 7 , F 1 3 . 4 , 2F 1 0 . 4 , I  16 , F 1 3. 4 , 2F1  C .4 , 19 ) 

45  CONTINUE 

LDA (  IT )  AJP3(5C)  ENI6(1B)  SLJ(IO) 

50  CONTINUE 
END 
END 
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APPENDIX  IV 
PROGRAM  FLOW  CHARTS 


SYMBOLS  USED  IN  THIS  APPENDIX 


START 


a  program  ^  subroutine  ,  box  ,  or 

page 


Decision  or  branch 


General  other  computer  operation 


Position  designation 


[uk-lmnI 

(A) 

(Q) 


132 


Logical  product  ITK  and  LMN 
Contents  o*f  the  "  A  register ss 
Contents  erf  the  n  Q.  register 

Statement  number  7S2. 
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PROGRAM  MASTER 


Boxes  I  B 


100 


Boxes  (o-\ 0 


101 


*  e>ox  12. 


102 


Boxes  \Z-!5 


103 


Boxes  14- 17 


104 


Bo*  17 


DO  MS 
Kr I ,KMAX 


ITEfAP=  NM6<) 


105 


Box  18 


106 


Bo  x  1 9 


107 


Box  20 


TEMP*  %  cold -tY1" 


i 


TEMprJ 
6(4  -  6, H 

12 

.  5  IS  ■  •  •  ^MASS  'TEW? 

l)=T£K?*C0LD(4-4>) 

108 


Box  21 


IN3=0 


Go  To  277 
(Box  -25) 


TLAST(M)  *  TTIME-TSL+TIME 
42 2(C)- 101*  PERCENT 
TEMP  =  B(7>  N\) 

LOCATE  =0 


Go  To  513 

(another  page) 


Clear  bit  *2  of  LbCM> 
ZTEMP  =  lO-ri  octal  positions 
TTtMP=6-4  octal  positions 


Clear  MASK2 
bits  of  LB(M) 


TT£MP= 

Set  bit 4*2 

TTEMP-Sf'5 

of  LBCM) 

Clear  bit  **3 
of  CBtlTEMP) 


Clear  MASKS  bits 
of-  L6(M) 

ITEM?  =6  it  43  of 
UZC^  •  C\eer 
bit  4S  of  UCM) 


Clear  MASK  I  blfc 
of  LfctTTEMP) 


Go  To  514- 

(anoihff  page) 


CATER.(W2)=  M 
W2.=WH+  I 


MISTAKEN)  = 
MrSTAKEC3)4  1 

& 


Clear  MASK! 
bits  of 
LBCM) 

GOO 


l 


MJ  -  M 
MOVE  =  I 
NOW  =0 
CALL  Tk/VAC  (mt) 
MOVE  -O 
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Go  To 

(a) 

Add  one  to 

k)e*t 

MISTAKE  (5) 

?a^e 

J 


Go  To  SI?) 
(next  pa<je) 


Put  TTEMP  in 

_ * _ 

*01 

lb(m) 

ITEMP=I-2'5 

JTEMP  =  M-3.,? 

I  =  NA  MV 

Put  JTEMP  in 

LB(NAV)^)-ELBCC> 

N  ext 


pacje 
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Ill 


e>o*  22 


112 


Boxes  ,23-25 


113 


Boxes  2C>~2S 


oo  *05 
J=|,IMMAX 


<o 

ENERGY =  ZCl.K) 

'ENERGY-  >-* - 

IEWEM.Y*  K  ■ 

i. 

>o 

305 

BCl-3,K>«B(l-3,K^ 
+3(4-4, K>  TSL 

_  \ 

[ 

Set  bit  *1 
MASK4  o-P 

and  cleat* 

LB(K) 

N\T=  LATERClT)  CALL  INVACCMJ) 
Upper  hal-f  c>4  IU  *  M 
Loner  hal-f  &  TO  =  M  J 
Write  Tape  6  ,  XU  >  6 


>o 


DO  241 

ITI  =KI2“  l 

261 

Go  To 

ir=  uiti 

Boy  2% 

t>0  79^  8 

DO  BOG 

T=1,INWAX 

K=I,XT| 

Put  KJFACE  m  octal  *1+  of  LVjLtti) 


799? 


NFACt=  0 


ITEMP=  octal  10-13  of  LB(fA') 

TTEhP=  octal  C-9  of  LB(M') 
Clear  MA$<*  bt U  of  Lfc(M> 


Clear  4B  bit 

C\ear  MASK\  blH 

o$  LBdTEMP'i 

""XTTEMP  >- — ► 

o£  LB(JTEMP) 

— KD 
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115 


Box  30 


116 


SUBROUTINE  SIT£(M) 


Box  31 


117 


SUBROUTINE  INVAUMJ) 


2>o  x  32 
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Box  33 


DO  2115 

y«\. 

>0 

JA=TAH 

1  =  1, 3 

- ►<IPACE(I)> 

TW=TA“l 

£0 


=o 


G> 


*0 

JC=I 

1 

k 

MIST  A  KE(2\  -  MI  VTAKE(  a-)  +- 1 

“ 

TPMAX=  S’ 


<2 


TPM  AX  =  ^ 

JPMAX=6 

1 

r  \ 

1 

r 

JPMAX=I3 

N00M(moVo 
H0ONi(l-\S)=O 
1=  1 


119 


120 


(A) 

wl  l\  cowtai  n 

2.451 

IL4-NLIME  +  TL 

2452 

IL-NLlbiE  +  TL 

2453 

XL  +  NPLANE  +  TL 

2.454 

XL-  NPLANJE  40-L„ 

2455 

IL  p 

2454 

XL- WLIKJE4TL+1. 

2457 

XL4KJPLAME4JL-H  _ 

245S 

XL+  NPLAKI.E4NLXNE 

245<* 

JL4WPLANE-NLXNE 

2440 

IL  -NPLANE+TL+I 

2441 

1L-NPLANE4NLXNE 

2442 

il-nplane-nlin/e 

2445 

XL~\  .  .  • 

2470 

IL4NLINE  43L  +  I 

n\ookj(x)=  itemp 
r  -  i  +  i 


TPMAX  =  I- 

-1 

Jp=l 


Noonctp)= 

MOONOO 


=o 


DO  2485 
K=I,TPMAX 

_ 3 

r 

X  = 

MOOWdxA 

121 


122 


123 


124 


Box  3^ 


125 


subroutine  inter. (mq,ii,s) 


Boxes  40-41 


126 


127 


Box-  -4E> 


128 


PROGRAM  SLAVE 


Box  I 


START 


REWIND  TAPE  8 


3 


/  f>£/\D  TAPE  8 

ix,  iy,i2,  abo-3),n/lia/£,  mplame.,maxlat,lmax, 

A,  SPHI,  FMASS,  ET/-/R£SH»  THERMAL,  EP&  0  LULL, 
ALPHA.*  SETA,  UTS,  MN/TS  ,  TTIME  ,TS  Li  Cl -30)  , 
A/TSCO-So'),  OBKJERGY  ,  TEV£R6,y  ,  MISTAKE n -  10), 
YEA/TRY  Cl-  10}  ,  ZEWTRY  0-103  ,  UA  ,  EMAX(  l- 30>, 
\^EFOb>MD>  LBC  l-LMAY),  BQ-7,  l-LW) 


Tupe  out  veminciev  to 
set  selective  switch  2.  up 


a 


RE  AO  1  CARO  (TAPE  2) 


0* 


DO  4-oo 
i=  1,4, 


PRINT  i  LINE 
(COMMEKJT  READ  IN) 


■Roo 


READ  XSY  KETC  l-^s) 


129 


130 


Box  2 


131 


Sox  3 


* 


0*i 


av\ot>hev 


pa^e 


132 


*  Next  pa^e 


133 


1445 


I 


Set  bit  *=4  of  Lg Cl) 


Go  "To 

^  *>  s' 

1 

Set  MASKS 0 

bits  of 

i-BO-} 

►  x= 

IL 

1 

r 

Set  bit*bl 
of  LB(I) 

&©x 

4 

NVAC=NVAC+MAXUT-  LMAX 

1456 

Add  one  to 

IMTERST 

134 


Box  4". 


(  PRINT:  IX,  AS(0, 
V  SPHI,  ,  ETMR 


135 


138 


137 


Box  G 


138 


Box  7 


=  0 


1  r 

-JL 


1517 


1521 


PRINT:  T^TEM, 
TEMP,  IC(l) 


i 


WUMBERS(  l-30^=IC(l-30) 

I L*  ~7 

CALL  HISTOfeo) 


1531 


PRINT:  I, IT, 
I,  IC(I+  36s) 


DO 

ISAI 

1  = 

1  >  30 

7=  1-7 
TEM  =  7-0.5 
TEMP*  7+0.5 


DO  1531 

x*  i , ao 

1 

IT =1-1 

139 


Box  8 


140 


»{VRINT:  NEUL? 


T  -  1-7 
TEM^  J-O.S 

TEMP  =  J+O.S 


(  PRINT:  LjE(l-7)iX> 
I  OCTALO-^  DEPTH, 
V  PERDTST 


IC(HOC»>=0 
IT  =  l 


Box  9 


DO  \S1%  L=t)LMAX 


Skip  to  next  \ 

pa<je  of  ouipuft  J 
papev'  J 


CALL  lKiITIAL(L) 
CALL  DEPTHPtL) 
CALL  ECOU tJT(L) 


PRINT:.  J, 
TEM^  TEMP, 
1C(I) . 


1574 


IL*0 


(  PRINT: 

\  ri  IT ,  I, 

Vic£l+sq)_ 


CALL  HlSToClo) 


IT  =  I-  | 
ICCT+B<rt 


PRINT:  I,  KETCiv 
KETft*i\lC(ic+«^, 
for  i*i>as 


NUMBERS^) 

*IC(5l-75i 


L 


Box 

CALL 

\0 

HIETO  (2$ 

141 


*  Mex"t  p3<^€ 


142 


143 


144 


SUBROUTINE  ROUND(L) 


Box  12 


145 


SUBROUTINE  INITLAL(L) 


box  13 


146 


SUBROUTINE  ECOUNT(L) 


Boxes  14  at*)  15 


1= 

1 

J 

L 

Ife40 

147 


SUBROUTINE  HISTO(NGROUPS) 


Box  I  £> 


148 


PROGRAM  RON 


149 


150 


151 


APPENDIX  V 


DESCRIPTION  OF  THE  FLOW  CHARTS 
The  explanation  given  below  is  to  be  used  with  APPENDIX  III 
(Program  Listings)  and  APPENDIX  IV  (Program  Flow  Charts) . 

PROGRAM  MASTER 

BOX  1 

This  box  writes  the  title  statement,  "THIS  IS  PROGRAM  MASTER" 
and  as  many  comment  cards  (80  spaces  per  card)  as  desired.  The  DO 
loop  is  written  for  five  comment  cards  at  present,  but  can  be  readily 
changed  if  desired. 

BOX  2 

Reads  input  parameters  for  the  run.  Formats  are  given  in  the 
listing  in  APPENDIX  III. 

BOX  3 

Changes  lattice  dimensions  (DC,  IY,  IZ)  into  floating  point  format, 
and  computes  the  physical  dimensions  of  the  lattice.  These  are  then 
permanently  found  in  AB  (1-3) .  It  must  be  kept  in  mind  that  IX  must 
be  ODD . 

BOX  4 

Computes  NLINE ,  NPLANE ,  MAXLAT ,  LMAX  (the  initial  LMAX)  from 
the  formulas  given  in  APPENDIX  II. 

BOX  5 

Computes  initial  lattice  positions  by  use  of  subroutine  SITE. 
Stores  the  site  number  and  the  fact  that  the  site  is  occupied  in  the  LB 
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array.  Stores  a  -10.0  in  each  location  of  the  TLAST  array.  This  is 
necessary  for  collisions  occurring  in  the  first  time  step#  otherwise 
these  collisions  would  be  ignored. 

BOX  6 

Writes  initial  lattice  positions  on  Tape  6  in  binary  for  eventual 
use  by  program  RON.  Also  sets  initial  or  zero  values  for  the  indicators 
shown. 

BOXES  7,8,  and  9 

Sets  initial  constant  values ,  sets  values  of  masks ,  and  reads  in 
the  input  parameters  for  the  appropriate  Born-Mayer  potential,  either 
potential  1,  2,  or  3. 

BOX  10 

This  box  includes  the  regeneration  process.  Jump  Switches  1  &  2 
are  raised,  leaving  Jump  Switch  3  down.  At  the  start  of  the  next  time 
step,  the  regeneration  information  will  be  written  on  Tape  7  in  binary, 
and  the  tape  will  be  rewound.  To  regenerate  the  run  the  program  must 
be  restarted,  with  all  three  jump  switches  in  the  raised  position.  After 
the  program  has  started  to  read  Tape  7  the  jump  switches  may  be  reset. 
When  the  computer  has  finished  reading  Tape  7 ,  the  tape  will  be  re¬ 
wound.  Further  regeneration  tapes  may  be  made  after  this  point  has 
been  reached.  The  main  DO  loop  over  the  time  steps  is  also  started  in 
this  box.  Immediately  after  the  start  of  a  time  step  the  current  values 
of  TENERGY,  N,  IENERGY,  and  ENERGY  are  written  on  Tapes  3  &  4. 

The  output  on  Tape  4  is  usually  transferred  to  the  typewriter  to  provide 
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a  visual  check  on  the  validity  of  the  run  as  it  progresses. 

BOX  11 

All  decisions  concerning  continuation  or  abandonment  of  the  run 
are  made  in  this  box.  A  jump  can  be  made  to  the  output  section,  the 
normal  output  dump  made,  and  a  return  to  this  box,  or  abandonment, 
accomplished  by  the  outside  control  of  the  jump  switches  o  The  run 
will  be  abandoned,  after  an  output  dump,  if  the  most  energetic  par¬ 
ticle  has  an  energy  less  than  THERMAL.  However,  the  run  will  not  be 
abandoned  under  this  condition  if  other  particles  are  to  be  shot  into  the 
lattice  o 

Additional  particles  will  be  shot  into  the  lattice,  if  the  proper  jump 
switch  is  set,  at  the  start  of  the  next  time  step  that  is  one  less  than  a 
multiple  of  LULL.  For  example,  if  LULL  is  3,  and  the  proper  switch  is 
set  in  time  step  6  (or  7) ,  an  additional  particle  will  enter  at  the  start 
of  time  step  8.  The  switch  must  be  left  in  the  set  position  until  after 
the  particle  has  been  read  into  the  program. 

The  input  co-ordinates  for  the  particle  are  given  as  -  1.0  for  X, 
and  the  actual  point  on  the  front  face  of  the  lattice  we  desire  to  hit. 

The  angles  ALPHA  &  BETA  are  given,  as  is  the  energy.  This  box  then 
calculates  the  proper  co-ordinates  and  velocity  components  to  achieve 
the  desired  impact  point.  Since  the  new  particle  is  the  most  energetic 
in  the  lattice,  the  Time  Step  Length  is  also  recalculated  at  this  time. 
BOX  12 

Here,  the  TSL  is  recalculated  if  the  number  of  the  time  step  is  an 
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even  multiple  of  JPB ,  and  the  first  30  values  are  placed  in  the  TSLI 
table  for  output  reference . 

BOX  13 

This  is  the  only  energy  check  made  during  the  run,  done  every  time 
step.  The  total  energy  (obtained  as  a  sum  of  the  individual  energies) 
is  compared  with  the  original  energy  (the  total  amount  shot  into  the 
lattice) ,  and  if  they  differ  by  more  than  1%,  the  MISTAKE  (1)  error 
indicator  is  increased  by  one, 

BOX  14 

The  value  of  ENERGY  is  recorded  in  the  EMAX  table  for  the  first  30 
time  steps,  the  total  energy,  TENERGY,  is  recalculated,  and  the  LB 
array  is  reset  and  now  shows  that  none  of  the  atoms  have  been  through 
the  current  time  step. 

BOX  15 

The  major  DO  loop  on  the  particles  is  started  here,  with  the  LB 
array  used  to  check  certain  conditions  on  each  particle.  If  the  atom: 

1)  Does  not  have  an  energy  greater  than  THERMAL,  or 

2)  Has  left  the  lattice ,  or 

3)  Has  completed  this  time  step,  or 

4)  Has  been  through  INVAC  and  has  not  had  a  collision  since , 
a  jump  is  made  to  the  end  of  the  DO  loop,  and  the  next  atom  is  con¬ 
sidered. 

Therefore,  only  those  atoms  are  considered  that  satisfy aU  the 
following  conditions: 
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1)  The  atom  is  a  bullet,  i.e.  has  more  energy  than  THERMAL, 

2)  The  atom  is  in  the  lattice , 

3)  The  atom  has  not  completed  this  time  step,  and 

4)  The  atom  has  not  been  through  INVAC  since  its  last  collision 
If  all  the  preceeding  conditions  are  satisfied,  the  initial  values  of  the 
counters  and  lists  shown  in  BOX  15,  APPENDIX  III  are  set  to  zero  or  other 
appropriate  values.  M,  the  fifth  Index  Register,  is  then  set  equal  to  L, 
the  fourth  Index  Register.  L,  as  such,  is  not  used  for  the  rest  of  the 
DO  loop,  until  this  pass  through  the  DO  loop  is  completed  and  a  new 

L  is  chosen . 

BOX  16 

All  atoms  within  a  mathematical  box  with  sides  3.02  (A)  (centered 
at  M)  are  placed  on  the  NUM  list  here.  Atoms  are  excluded  if  they  have 
left  the  lattice  or  have  been  through  the  time  step.  If  the  list  contains 
only  one  atom  (M  is  placed  first  on  the  list) ,  then  a  jump  is  made  to  BOX 
25  where  another  value  of  M  is  chosen  from  the  LAST  list.  With  other 
atoms  in  the  box,  indicators  IN3  and  NFACE  are  set  to  zero. 

BOX  17 

For  every  atom  on  the  NUM  list,  the  following  parameters  are  com¬ 
puted  with  respect  to  atom  M: 

1)  The  undeviated  distance  of  closest  approach  between  centers , 

called  DSTANCE  (  ) , 

2)  The  time  of  closest  approach  between  centers,  also  un¬ 
deviated,  called  TMIN  (  ), 
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3)  The  time  when  the  two  atoms  are  separated  by  a  distance  of 
SPHI,  called  T  (  ).  This  is  the  actual  time  a  collision 
would  occur. 

The  distance  is  given  in  Angstroms ,  with  the  time  in  jiffys  ,  relative 
to  the  start  of  the  current  time  step. 

Various  conditions  are  now  imposed  on  these  quantities  to  insure  the 
validity  of  these  presumed  collisions .  An  atom  must  pass  all  the  con¬ 
ditions  to  insure  that  a  collision  will  take  place.  The  conditions  are: 

1)  &V)  >  1.0  x  10"6, 

2)  DSTANCE  (  )  SPHI , 

3)  T  (  )  lies  between  +  TSL  of  the  start  of  this  time  step, 

4)  T  (  )  is  less  than  0.99999  TSL.  This  eliminates  collisions 
that  end  too  far  into  the  next  time  step . 

5)  The  TLAST  for  both  atoms  must  be  less  than  the  absolute 
time  of  this  collision. 

If  all  of  the  above  conditions  are  satisfied,  IN3  is  increased  by 
one.  If  one  or  more  of  the  conditions  is  not  satisfied,  T  (  )  is  set 

equal  to  a  large  positive  number  (100  TSL  serves  this  purpose  in  the 
program) . 

It  should  be  noted  that  DSTANCE  (4) ,  TMIN  (4) ,  and  T  (4)  are  the 
parameters  associated  with  the  atom  that  is  number  4  on  the  NUM  list, 
i.e.  NUM  (4),  and  so  on. 

BOX  18 

The  minimum  time  on  the  T  (  )  list  is  found.  All  those  on  the  NUM 
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list  whose  corresponding  T's  are  equal  to  the  minimum  time,  or  within 
(CUTOFF) (TSL)  of  the  minimum  time  are  than  placed  on  the  KHIT  list. 

Atom  M  is  placed  first  on  the  KHIT  list.  At  this  point,  if  INPUT  is  zero, 
the  KHIT  list  is  duplicated  on  the  LAST  list.  This  occurs  only  when  M 
equals  L.  INPUT  is  increased  and  the  next  box  is  entered  if  there  is 
more  than  one  atom  on  the  KHIT  list,  i.e.  another  atom  besides  M.  If 
there  are  no  other  atoms  on  the  KHIT  list,  a  jump  is  made  to  BOX  25, 
and  M  is  set  equal  to  the  next  atom  on  the  LAST  list. 

BOX  19 

This  box  and  the  next  are  concerned  with  the  two  body  interaction 
itself,  and  the  subsequent  scaling  methods  for  energy  and  velocity. 

These  scaling  methods  are  used  whenever  M  collides  with  more  than 
one  atom.  The  pre -interaction  velocities  and  energies  are  stored  in 
the  SAVE  array.  The  COLD  list  is  set  to  zero,  and  J  is  set  to  KHIT  (JA) . 
JA  is  the  index  used  for  the  KHIT  list  in  this  box,  M  then  hits  the  atoms 
in  order  of  their  placement  on  the  KHIT  list.  JA  is  initially  two,  and 
increases  up  to  and  including  NMAX. 

With  J  selected,  LB  (J)  is  checked,  and  if  J  is  a  member  of  an  inter¬ 
stitial  pair,  NPAIR  is  set  to  3.  The  indicator  NPAIR  is  used  in  BOX  24 
to  insure  that  the  appropriate  action  is  taken  with  respect  to  J's  inter¬ 
stitial  partner  after  the  collision  is  over. 

Since  the  atoms  on  the  NUM  and  KHIT  lists  do  not  have  the  same 
placement,  a  search  is  now  made  to  find  J  on  the  NUM  list.  When  this 
is  done,  we  have  the  collision  time,  T,  and  the  collision  parameter,  PAR. 
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The  two  atoms ,  M  &  J ,  are  then  moved  ahead  in  time  to  time  T .  If  M 
has  already  had  one  or  more  collisions,  it  must  be  moved  back  in  time 
to  the  start  of  this  time  step,  and  then  ahead  to  the  start  of  this  collision. 

The  interaction  subroutine  is  then  called.  Outputs  are  contained  in 
the  DEV  array  as  the  new  positions,  velocities,  and  energies  of  M  &  J. 
TLAST  is  now  calculated  for  J,  and  J's  co-ordinates  and  velocities 

are  set  to  those  of  the  DEV  array.  Notice  that  the  TLAST  for  J  does  not 

% 

include  any  portion  of  the  time  necessary  for  the  collision.  If  a  portion 
of  this  time  were  included,  significant  collisions  would  be  eliminated. 

The  COLD  array  (1-6)  contains  the  sum  of  the  final  position  and  velocity 
co-ordinates  for  M  as  given  by  the  subroutine  for  all  of  M's  collisions. 
The  total  energy  lost  by  M  (ELOST)  then  equals  the  energy  gained  by 
this  J  plus  that  gained  by  the  other  J’s  on  the  KHIT  list.  NHIT  shows 
the  actual  number  of  atoms  that  have  collided  with  M . 

BOX  20 

In  this  box,  the  scaling,  by  various  means, of  the  velocities  and 
energies  of  the  atoms  on  the  KHIT  list  is  accomplished. 

The  temporary  quantity  T2  is  set  equal  to  the  previous  energy  of  M 
minus  the  amount  lost  by  M  in  all  its  collisions . 

If  T2  is  zero  or  negative*,  the  energy  and  velocities  of  M  are  set 
to  zero.  With  T2  negative,  the  energies  and  velocities  of  all  others 
on  the  KHIT  list  must  be  scaled  to  conserve  energy.  All  atoms  are 
scaled  in  the  same  proportion,  the  absolute  direction  of  their  velocity 
vectors  is  preserved ,  their  energies  being  scaled  by  the  square  of 
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the  scale  factor  used  on  the  velocity  components  . 


In  all  cases,  the  spatial  co-ordinates  of  M  after  the  collision 
must  be  found.  These  are  set  equal  to  the  numerical  average  of  M’s 
positions  as  computed  by  the  subroutine.  This  is  accomplished  by 
simply  dividing  COLD  (1-3)  by  NHIT . 

For  a  positive  T2,  there  are  two  cases ,  the  energy  of  M  before  the 
collisions  above  ETHRESH  and  below  ETHRESH. 

When  the  original  Ej^  ETHRESH  and  there  was  only  one  collision, 
the  co-ordinates  of  M  are  set  to  the  DEV  results,  and  a  jump  is  made 
to  the  next  box.  If  there  was  more  than  one  collision,  Em  is  set  equal 
to  T2,  and  the  velocities  are  scaled  so  that  the  components  will  be  in 
the  same  proportion  as  the  sum  of  the  individual  resultant  velocities 
as  computed  by  the  subroutine.  This  essentially  gives  M  a  resultant 
direction  equal  to  the  vector  sum  of  the  directions  M  would  have  as  a 
result  of  its  collisions  with  the  rest  of  the  KHIT  list. 

When  the  original  Em«^  ETHRESH,  we  assume  that  M  has  been 

oscillating  about  a  site,  and  its  site  number  is  obtained  from  LB  (M) . 

Em  is  then  set  equal  to  T2,  and  the  velocities  are  scaled  so  that  the 

resultant  velocity  vector  will  point  toward  M’s  oscillatory  center.  The 

scaling  factor  is  (V^/  (AX^  +AY^  +  AZ^))^?  This  is  multiplied  by 

&X,  Ay •  and  AZ  to  obtain  the  correct  components.  This  scaling  auto- 
2 

matically  keeps  V  in  proper  proportion  to  Em. 

BOX  21 

TLAST  for  M  is  calculated  at  the  start  of  this  box.  M's  TLAST 
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could  not  be  calculated  before  this  point  since  M  had  not  yet  parti¬ 
cipated  in  all  its  possible  collisions. 

This  box  then  insures  that  all  atoms  with  energies  less  than 
ETHRESH  have  site  numbers  stored  in  LB,  that  they  are  placed  on  the 
LATER  list  if  their  energies  are  less  than  THERMAL,  or  that  they  are  left 
alone  if  in  between  the  two  thresholds  . 

If  an  atom  has  E^  ETHRESH,  its  site  number  must  be  removed  from 
LB  (M) ,  and  the  bit  signitifying  an  occupied  site  must  be  removed  from 
LB  (site) ,  unless  the  atom  was  once  part  of  an  interstitial  pair.  If  it 
was  part  of  an  interstitial,  then  the  partner  section  of  LB  must  be  re¬ 
moved  from  both  LB's. 

If  an  atom  should  occupy  a  site  and  doesn’t,  then  a  site  must  be 
found.  This  is  accomplished  by  going  part  way  through  INVAC ,  INVAC 
is  exited  by  use  of  the  MOVE  indicator.  If  a  site  is  calculated  for  M, 
and  it  is  occupied,  the  information  pertinent  to  an  interstitial  pair  is 
recorded  in  LB.  If  the  site  is  occupied  by  an  interstitial,  M  is  placed 
on  the  site  anyway,  and  the  fact  that  a  "triple  interstitial"  was  formed 
is  recorded  by  increasing  MISTAKE  (5).  Em  is  inspected  at  this  point, 
and  MISTAKE  (3)  is  increased  if  Em  is  negative.  The  symbolic  language 
of  this  box  may  appear  complicated,  but  comprehension  should  be 
facilitated  by  the  flow  chart. 

BOX  22 

All  other  atoms  on  the  KHIT  list  are  now  subjected  to  the  same  pro¬ 
cess  M  underwent  in  the  previous  box,  except  that  here  a  jump  is  made 
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to  Statement  600  to  enter  INVAC .  To  preserve  the  values  of  I ,  J,  and  M, 
and  also  make  M  equal  J,  the  values  of  I,  J,  and  M  are  stored  in  MAD, 
NAVY,  and  LOCATE,  respectively.  These  two  boxes  might,  at  some 
future  stage  of  development  be  combined  into  one  box  covering  the  entire 
KHIT  list.  However,  when  this  section  was  developed,  the  treatment  of 
M  and  the  others  on  the  KHIT  list  was  sufficiently  different  so  that  this 
combination  would  have  been  prohibitively  difficult.  Most  of  this  dif¬ 
ficulty  has  disappeared  as  the  program  developed. 

BOX  23 

All  atoms  on  the  KHIT  list  are  now  backed  up,  time-wise  and  space- 
wise,  to  the  beginning  of  the  current  time  step-  Each  atom  must  be 
backed  up  using  its  particular  T  (  ) ,  since  they  were  not  all  hit  at 

exactly  the  same  time.  A  search  through  the  NUM  list  accomplishes 
this,  giving  the  correct  value  on  the  T  list  for  the  computation. 

BOX  24 

If  one  of  the  J’s  on  the  KHIT  list  was  a  member  of  an  interstitial 
pair,  then  NPAIR  is  three,  and  a  jump  is  made  to  Statement  236.  The 
relations  between  atom  M  (the  incoming  atom) ,  atom  J  (the  member  of 
the  pair  struck  by  M) ,  and  the  partner ,  atom  LL  (the  other  member  of  the 
interstitial  pair) ,  are  now  considered. 

If  J  has  an  energy  between  the  two  thresholds  (vibrating  about  its 
site),  the  information  in  LB  about  it  must  be  retained.  If  the  partner, 

LL,  has  been  hit,  and  has  an  energy  greater  than  ETHRESH,  the  MASK1 
portion  of  both  LB’s  must  be  erased.  If  both  have  an  energy^  ETHRESH, 
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nothing  is  done  to  either  LB.  If,  however,  J  has  an  energy  greater  than 
ETHRESH  and  atom  LL  has  not  been  hit,  we  must  erase  from  LB  (J)  the 
site  and  partner  number,  erase  from  LB  (LL)  LL's  partner  number  (which 
is  J) ,  and  also  move  atom  LL  back  onto  the  site.  It  would  be  a  definite 
error  to  leave  the  partner,  LL,  on  its  previous  split  position,  so  it  seems 
more  realistic  to  move  it  back  onto  the  site. 

BOX  25 

The  interactions ,  scaling ,  and  correction  of  LB  have  now  been  com¬ 
pleted  for  this  round  of  collisions.  A  jump  is  now  made  to  Statement  200 
in  BOX  16.  If  M  has  no  more  collisions  this  time  step,  a  jump  back  to 
this  point  will  be  made.  If  M  does  have  further  collisions,  they  are  cal¬ 
culated  and  carried  out. 

Upon  entering  this  box,  if  IN3  is  zero,  then  no  collisions  were  made 
during  the  last  pass  through,  and  this  M  has  completed  the  time  step. 
The  next  atom  on  the  LAST  list,  the  INMIN  th  one, is  made  atom  M. 

INMAX  is  the  limit  of  this  list.  INMIN  is  increased  by  one,  and  a  jump 
is  made  to  Statement  200.  When  the  LAST  list  has  been  exhausted,  i.e. 
each  one  has  been  through  the  process  until  it  will  have  no  more  colli¬ 
sions  this  time  step,  an  exit  to  the  next  box  is  made. 

Any  atoms  on  the  LAST  list  that  have  been  put  on  the  LATER  list  to 
go  through  INVAC  are  not  sent  back  through  the  collision  process  since 
they  have  MASK6  set. 

The  name  of  the  atom  under  primary  consideration  was  changed  from 
L  to  M  in  BOX  15  to  prevent  the  above  process  from  destroying  the  value 
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of  L,  since  L  is  the  index  on  the  main  DO  loop  over  the  atoms . 

BOX  26 

In  a  DO  loop  over  the  LAST  list,  only  those  atoms  not  on  the  LATER 
list  are  selected,  and  then  moved  to  the  end  of  the  time  step.  This  fact 
is  then  recorded  in  LB.  The  atom's  energy  is  then  compared  with 
ENERGY,  the  largest  one  becomes  the  new  ENERGY,  and  INERGY  is 
changed  if  needed. 

BOX  27 

If  NZ  =  1 ,  the  LATER  list  is  empty,  and  a  jump  is  made  around  this 
box.  If  not,  IT1  is  set  equal  to  (NZ  -  1) ,  as  the  upper  limit  of  the  DO 
loop  to  follow.  INVAC  is  then  called  successively  for  all  of  the  LATER 
list,  writing  on  Tape  6  in  binary  the  particle  number  and  time  step, 
both  contained  in  IU,  and  the  X,  Y,  and  Z  co-ordinates. 

BOX  28 

In  this  DO  loop  on  J,  only  those  atoms  on  the  LAST  list  that  were 
not  on  the  LATER  list  (and  have  consequently  been  through  INVAC)  are 
considered . 

All  three  co-ordinates  are  compared  to  the  lattice  dimensions .  If 
an  atom  is  outside  the  physical  dimensions  of  the  lattice,  a  value  of 
NFACE  is  calculated.  NFACE  and  the  number  of  the  time  step  of  exit 
are  then  recorded  in  LB  (M) ,  the  site  number  and  possible  interstitial 
partner  are  erased  from  LB  (M) ,  and  M's  number  is  erased  from  any 
interstitial  partner's  LB.  The  memory  of  occupancy  is  erased  from  the 
site  if  M  is  not  an  interstitial  member,  and  the  next  atom  on  the  LAST 
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list  Is  considered. 


This  is  the  end  of  the  major  DO  loop  on  L,  the  next  L  is  chosen, 
or  an  exit  is  made  to  the  next  box  if  the  DO  loop  on  L  is  completed, 

BOX  29 

The  same  procedure  is  followed  here  as  was  followed  in  BOXES  26 
and  28  above ,  except  that  atoms: 

1)  That  have  been  through  the  time  step , 

2)  Having  an  energy  less  than  THERMAL,  or, 

3)  That  are  outside  the  lattice 
are  not  considered. 

This  is  the  end  of  the  major  DO  loop  on  N,  the  time  step  is  completed, 
and  either  a  new  time  step  is  started ,  or  an  exit  to  the  next  box  and  the 
output  section  is  made. 

BOX  30 

Writes  selected  lattice  and  program  parameters ,  and  a  general  lattice 
dump  on  Tape  8  in  binary  for  eventual  use  by  the  SLAVE  program.  Various 
BCD  output  is  then  written  on  Tape  3 ,  and  if  the  DO  loop  on  N  (the  time 
steps)  was  completed,  a  (-1)  is  written  in  binary  on  Tape  6  and  the  pro¬ 
gram  ends.  If,  however,  the  DO  loop  on  N  was  not  completed,  a  jump  is 
made  back  to  BOX  11.  There,  a  jump  is  made  back  to  this  box  (with  the 
resultant  (-1)  on  Tape  6,  etc.)  if  Jump  Switch  3  is  not  set.  In  the  normal 
operational  mode,  the  DO  loop  is  not  completed,  and  the  program  itself 
terminates  computation  by  a  jump  from  BOX  11  to  the  start  of  this  box, 
the  output  is  written,  the  jumps  back  to  BOX  11  and  back  are  made,  and 
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the  program  ends . 

A  jump  around  the  BCD  section  of  the  output  may  be  made  by  setting 
Jump  Switches  1  &  3.  To  facilitate  resetting  of  the  jump  switches,  Stop 
Switch  1  may  be  set.  If  this  is  done,  the  program  will  stop  after  writing 
Tape  8 ,  and  again  before  jumping  back  to  BOX  1 1 . 

Tapes  3  &  4,  then,  are  used  for  BCD  output,  and  Tapes  6,7,  and  8 
are  used  for  binary  output  over  the  course  of  a  run  (including  regeneration) . 
Since  Tape  4  output  is  normally  switched  to  the  typewriter,  only  four  out¬ 
put  tapes  are  needed  in  the  course  of  a  run.  The  three  binary  tapes 
(6,7,8)  will  be  rewound  at  the  end  of  the  run,  and  may  be  used  without 
further  action .  An  END  OF  FILE  mark  must  be  placed  on  the  end  of  Tape  3 
before  printing. 

SUBROUTINE  SITE  (M) 

BOX  31 

This  short  subroutine  calculates  the  fixed  point  co-ordinates  of 
the  lattice  site  given  as  ah  input  parameter.  The  co-ordinates  are  found 
in  IB  (1-3)  after  exit  from  the  subroutine. 

SUBROUTINE  INVAC  (MJ) 

This  subroutine,  given  the  number  of  an  atom,  either  places  the  atom 
on  a  vacant  lattice  site  near  it,  or  forms  an  interstitial  pair  with  one  of 
its  nearest  neighbors . 

BOX  32 

The  initial  indicator  IN 2  is  set  to  zero,  and  the  input  atom  is 

v 

changed  to  number  MM.  The  site  number  occupied  by  MM,  and  the  fact 
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that  MM  is  on  the  LATER  list  are  now  cleared  from  LB  (MM),  and  the 
fact  that  MM  has  gone  through  INVAC  is  set  in  LB  (MM) . 

Atom  MM's  co-ordinates  are  now  inspected  to  see  whether  or  not  MM 
is  within  one  unit  of  a  face  of  the  lattice.  If  so,  the  number  of  this  face 
is  stored  in  IFACE  (1-3) ,  and  the  indicator  NIF  is  increased.  The  index 
of  IFACE  is  1 ,  2,  or  3,  depending  on  whether  the  face  is  in  the  X,  Y,  or 
Z  direction,  respectively. 

The  co-ordinates  of  MM  are  rounded  off,  and  the  fixed  point  co¬ 
ordinates  of  the  nearest  site  (or  non-site)  are  found. 

BOX  33 

We  must  now  determine  which  of  the  two  general  cases  we  have,  site 
or  non-site,  and  if  IFACE  (1-3)  is  non-zero,  which  of  the  12  special 
cases  we  are  to  consider.  If  the  sum  of  the  three  fixed  point  co-ordinates 
is  even,  then  they  designate  a  lattice  site,  and  JM  is  set  to  1 .  If  the 
sum  is  odd,  then  a  non-site  is  designated,  and  JM  is  set  to  zero.  The 
sum  of  the  Y  and  Z  co-ordinates  is  then  checked.  If  odd,  JL  =  0,  and  if 
even,  JL  =  -1.  If  the  indicator  MOVE  is  non-zero,  a  RETURN  to  the  main 
program  is  executed. 

The  formula  for  IL  (the  site  number)  results  in  an  X  co-ordinate 
one  less  or  one  more  than  an  actual  site  for  the  non-site  case.  The  sum 
of  JL  and  JM  determine  which  case.  If  the  sum  is  -1,  then  one  is  added 
to  IL. 

At  this  point,  IL  is  always  either: 

1)  The  site  of  the  rounded  off  position  of  MM,  or 
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2)  The  site  with  X  co-ordinate  one  larger  than  MM's  rounded-off 
position. 

The  possible  special  case  is  now  determined.  If  IFACE  (1-3)  ^  0,  then 
JB  is  set  equal  to  the  face  number  with  JD  =  1,  If  more  than  one  IFACE 
is  non-zero,  MISTAKE  (2)  is  increased,  since  this  means  a  position  along 
an  edge  or  comer  of  the  lattice . 

In  either  case,  if  IFACE  (1-3)  9^0,  then  JA  is  non-zero,  and  JC  is 
set  equal  to  1.  If  IFACE  (1-3)  is  zero,  then  JC  =  JM  (either  1  or  0) . 

JB  is  set  equal  to  the  first  non-zero  IFACE  or  to  zero  if  IFACE  (1-3) 
is  zero.  JD  is  then  set  to  1  if  a  non-zero  IFACE  exists,  and  to  zero 
otherwise . 

The  case  number  ( 1  —  1 4)  is  now  established  by  the  formula: 

JW  =  (JM*6*JD)  +  1  +  JB  +  JC  .  Cases  1  &  2  are  the  general  cases. 

Cases  3-8  are  special  forms  of  case  1,  with  cases  9-14  as  special 
forms  of  case  2.  These  special  cases  simply  eliminate  some  of  the 
nearest  neighbors  found  in  the  general  case  from  consideration. 

The  lists  MOON  and  NOON  are  set  to  zero ,  and  the  next  box  is 
entered . 

BOX  34 

The  first  section  (down  to  Statement  2451)  is  a  section  of  computed 
GO  TO  statements  to  insure  the  correct  choice  of  the  nearest  neighbor 
formulas.  The  second  section  places  the  site  numbers  of  the  nearest 
neighbor  positions  on  the  MOON  list,  provided  their  numbers  lie  between 
one  and  MAXLAT.  I  is  the  index  used  for  the  MOON  list,  and  is  always 
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one  greater  than  the  number  of  sites  on  the  list, 

BOX  35 

All  of  the  nearest  neighbor  numbers  found  may  not  have  been  placed 
on  the  MOON  list  (possible  negative  numbers,  or  numbers  greater  than 
MAXLAT) ,  JPMAX,  therefore,  is  reset  to  (i-l).  If  the  indicator  MOVE  is 
non-zero,  a  RETURN  to  the  main  program  is  executed. 

The  LB  for  each  atom  whose  site  is  on  the  MOON  list  is  now  in¬ 
spected.  If  the  LB  shows  that  the  site  is  not  occupied  (see  APPENDIX  I) , 
the  site  number  is  placed  on  the  NOON  list.  JP  is  the  index  used  for  the 
NOON  list. 

The  MOON  list,  then,  is  a  list  of  the  sites  of  the  nearest  neighbors, 
and  the  NOON  list  is  a  list  of  those  nearest  neighbor  sites  that  are  vacant 
sites,  with  JK  the  number  of  vacancies  on  the  NOON  list. 

BOX  36 

If  JK  =  1,  then  only  one  vacancy  exists,  and  atom  MM  is  placed 
there,  the  site  number,  IL,  is  set  in  LB  (MM),  and  the  fact  that  the  site 
is  now  occupied  is  stored  in  LB  (IL) .  If  IN2  is  zero,  a  RETURN  is  then 
executed,  otherwise  a  jump  to  BOX  39  (Statement  2800)  is  executed. 

If  JK  is  zero,  there  are  no  vacancies  and  a  prospective  interstitial 
partner  must  be  found.  If  IN2  is  zero,  a  jump  to  BOX  38  is  made  to 
accomplish  this.  If  IN2  is  non-zero,  then  MM  is  the  prospective  inter¬ 
stitial  partner  of  the  original  atom  MJ ,  and  a  jump  is  made  to  BOX  39  to 
form  the  interstitial. 
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If  JK  >  1 ,  then  more  than  one  vacancy  exists .  A  jump  is  made  to 
BOX  37,  where  one  of  the  vacancies  is  chose,  a  return  to  this  box  is  then 
made,  and  the  procedure  given  for  JK  =  1  is  followed. 

BOX  37 

The  (distance) ^  between  atom  MM  and  the  vacant  sites  on  the  NOON 

2 

list  is  computed  and  stored  in  the  DSTANCE  array.  The  minimum  (distance) 
is  then  found,  and  the  DSTANCE  list  is  inspected  to  find  the  number  of 
sites  that  are  at  this  minimum  distance  from  atom  MM.  If  there  is  only 
one,  IL  is  set  equal  to  that  site,  and  a  jump  is  made  to  Statement  2492  in 
BOX  36 .  If  there  is  more  than  one ,  MM  is  placed  on  its  own  site ,  if  its 
own  site  is  one  of  those  on  the  list.  If  not,  then  MM  is  placed  on  the 
first  or  second  site  on  the  NOON  list,  depending  on  whether  the  number 
of  the  time  step  is  odd  or  even,  respectively. 

The  number  of  the  time  step  is  used  here  as  a  convenient  random 
number  generator,  since  the  decision  is  made  on  oddness  or  evenness, 
rather  than  the  actual  value  of  the  number  itself. 

BOX  38 

This  box  is  entered  only  from  BOX  36,  and  then  on  the  condition  that 
IN2  is  zero  and  that  MM  (MJ)  has  found  no  vacancies.  A  prospective 
interstitial  partner  must  be  chosen  from  the  MOON  list. 

The  time  of  atom  MM's  (MJ's)  closest  approach  to  each  of  the  sites 
listed  on  the  MOON  list  is  found  and  stored  in  the  TMIN  list.  Any  sites 
occupied  by  interstitial  pairs  are  not  considered. 
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The  minimum  time  of  closest  approach  is  then  found,  TMIN  (KMIN) , 
which  locates  the  prospective  interstitial  site.  MAX  is  derived  from 
A  X,  AY,  or  AZ,  whichever  is  greatest  (Z\  X  =  X  co-ordinate  of  MM 
minus  X  co-ordinate  of  the  site  chosen) ,  and  this  becomes  the  prospective 
axis  of  the  interstitial  pair.  The  atom  occupying  this  site  is  then  found. 
Since  it  is  possible  for  this  atom  to  be  pushed  into  any  vacancy  near  it, 
and  MM  to  fall  into  the  vacancy  thus  created,  this  atom's  nearest  neigh¬ 
bors  must  be  inspected.  This  is  done  by  setting  MM  equal  to  the  number 
of  this  atom,  setting  IN2  to  1,  and  jumping  to  Statement  2361  in  BOX  32. 
BOX  39 

If  the  new  MM  (the  original  atom  is  still  designated  by  MJ)  does  find 
a  vacancy,  it  is  placed  there  in  BOX  36,  a  jump  is  made  to  this  box,  MM's 
now  vacant  site  is  named  as  the  prospective  site  for  MJ,  a  return  to 
BOX  36  is  made,  and  MJ  is  placed  on  its  prospective  site.  A  RETURN  to 
the  main  program  is  then  executed  in  BOX  36. 

If,  however,  the  new  MM  did  not  find  a  vacant  site,  then  a  jump  is 
made  to  this  box  (Statement  2720),  and  the  interstitial  is  actually  formed. 
The  site  of  the  interstitial  is  entered  in  each  LB,  together  with  the  other 
atom’s  number.  The  fact  that  the  site  is  occupied  is  then  entered  in  the 
LB  of  the  site,  and  a  RETURN  is  executed. 

If  the  new  MM  equals  MJ,  then  the  atom  has  formed  an  interstitial 
with  itself,  and  the  error  counter  MISTAKE  (4)  is  increased.  The  inter¬ 
stitial  is  formed  so  that  MJ  is  on  the  same  side  as  it  was  before  forma¬ 
tion  of  the  interstitial  pair. 
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SUBROUTINE  INTER  (M Q,  II,  S) 

This  subroutine,  given  the  atom  under  consideration  (MQ) ,  the 
target  atom  (II) ,  and  the  impact  parameter  (S ,  in  meters)  calculates  the 
positions  of  the  atoms  at  the  end  of  the  collision,  their  velocities,  and 
energies  and  stores  them  in  the  DEV  array » 

BOX  40 

Transforms  atom  MQ's  position  co-ordinates,  velocity  components, 
and  energy  and  the  angle  between  the  radius  vector  from  atom  II  to 
atom  MQ  and  MQ's  velocity  into  the  target  (atom  II)  frame  of  reference. 
Sets  CPA  =  SPHI  as  the  initial  value  for  the  CPA  iteration. 

BOX  41 

Computes  CPA  by  an  iterative  process.  The  iteration  is  terminated 
when  two  successive  values  differ  by  less  than  10“^  % ,  or  when  any 
value  is  greater  than  the  preceeding  value. 

BOX  42 

Computes  the  time  required  for  the  interaction  and  the  total  angle 
of  deviation.  The  integrals  are  solved  by  numerical  integration  using  a 
Gauss -Legendre  four  point  quadrature. 

BOX  43 

Computes  the  composite  factors  and  trigonometric  functions  that  are 
used  more  than  once  in  the  transformation  to  the  laboratory  reference 
frame . 

BOX  44 

Computes  direction  cosines,  relative  to  the  laboratory  frame  of 
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reference,  of  a  line  in  the  plane  of  the  interaction  that  is  perpendicular 
to  the  line  between  the  two  atoms  at  the  beginning  of  the  interaction. 

BOX  45 

Computes  the  position  co-ordinates,  velocity  components,  and 
energies  of  both  atoms  and  stores  them  in  the  DEV  array.  See  APPENDIX  I, 
DEV  for  the  actual  location. 

PROGRAM  SLAVE 

BOX  1 

Reads  binary  input  from  Tape  8.  "SET  SLJ  SWITCH  2  UP  IF  NEEDED" 
is  written  on  Tape  3  (normally  transferred  to  the  typewriter) .  As  many 
comment  cards  (80  spaces  per  card)  as  desired  (6  at  present)  are  read 
from  Tape  2  and  written  on  the  output  tape .  Reads  input  parameters  for 
the  run  from  Tape  2.  Writes  on  the  output  tape  "SELECTIVE  JUMP  SWITCH 
NUMBER  TWO  WAS  UP"  if  Jump  Switch  2  is  used  for  the  run. 

IHISTO  (1-4)  and  IHISTO  (105)  are  defined,  the  IHISTO  array  then 
serves  as  the  format  specification  for  the  histographs.  The  constants 
I2E15,  I2E24,  F2E24,  and  I2E30  are  set,  and  the  masks  for  use  with  the 
LB  array  are  defined . 

BOX  2 

If  NTS  (called  N  in  PROGRAM  MASTER,  the  number  of  time  steps 
completed)  is  zero,  then  NTS  is  set  equal  to  MNTS;  other  constants  and 
indicators  are  then  set.  The  times  in  the  TSLI  array  and  the  time  TSL 
are  converted  from  jiffy s  to  seconds.  Various  portions  of  the  LB  array 
are  then  cleared,  depending  on  the  settings  of  the  jump  switches. 
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BOX  3 


In  this  box  all  atoms  are  placed  in  one  of  the  following  classifica¬ 
tions: 

1)  Atoms  occupying  their  original  sites  , 

2)  Atoms  occupying  sites  other  than  their  original  sites  , 

3)  Atoms  that  are  parts  of  interstitial  pairs , 

4)  Atoms  whose  energy  is  greater  than  ETHRESH,  and 

5)  Atoms  outside  the  lattice. 

Program  MASTER  has  also  calculated  the  face  of  exit  for  those  atoms  that 
have  left  the  lattice . 

The  number  of  atoms  and  the  total  energy  of  those  atoms  is  calculated 
for  each  classification,  as  is  the  number  of  vacant  lattice  sites. 

BOX  4 

The  total  number  of  atoms  and  their  total  energy  is  calculated. 

This  serves  as  an  additional  total  energy  check  at  the  end  of  the  compu¬ 
tational  run.  These  numbers  and  various  other  data  (see  APPENDIX  III) 
are  written  on  the  output  tape.  The  number  of  atoms  and  their  total 
energy  is  written  on  the  output  tape  for  each  classification  listed  in 
BOX  3. 

BOX  5 

The  number  of  atoms  occupying  their  original  sites  is  written  on  the 
output  tape.  An  array  of  these  atoms  is  then  written,  in  order,  on  the  out¬ 
put  tape.  Atoms  not  on  their  original  sites  are  represented  by  zeros  in 
this  array. 
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BOX  6 


The  total  number  of  vacant  lattice  sites  and  the  numbers  of  these 
sites  are  written  on  the  output  tape. 

BOX  7 

The  total  number  of  replacements  is  written  on  the  output  tape. 

The  number  of  the  atom,  the  number  of  its  final  site,  its  penetration, 
and  its  radial  penetration  are  written  on  the  output  tape  for  each  of 
the  replacement  atoms.  Histographs  are  then  compiled  and  written  on  the 
output  tape  for  the  penetration  and  the  radial  penetration  distributions 
of  the  replacement  atoms . 

BOX  8 

The  total  number  of  interstitial  pairs  is  written  on  the  output  tape. 
The  number  of  each  atom  in  a  pair,  its  partner,  the  site  they  occupy,  the 
penetration  of  each  atom  in  the  pair,  and  the  radial  penetration  of  each 
atom  in  the  pair  is  written  on  the  output  tape.  This  is  written  twice  for 
each  interstitial  pair  since  there  are  two  atoms  in  the  pair.  Histographs 
are  compiled  and  written  on  the  output  tape  for  the  penetration  and  the 
radial  penetration  distributions  of  all  interstitial  pair  members. 

If  Jump  Switch  2  was  set  for  the  run,  the  interstitials  are  determined 
by  inspecting  the  LB  array  given  by  program  MASTER.  If,  however,  Jump 
Switch  2  was  not  set,  the  program  disregards  the  LB  array  and  calculates 
the  interstitial  pairs  in  the  lattice . 

BOX  9 

The  total  number  of  atoms  with  energy  greater  than  ETHRESH 
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(knock-ons)  is  written  on  the  output  tape .  The  number  of  each  atom ,  its 
co-ordinates,  velocity  components,  energy,  its  LB  (in  octal),  its  pene¬ 
tration  ,  and  its  radial  penetration  are  written  on  the  output  tape  for  each 
knock-on  atom.  Histographs  are  then  compiled  and  written  on  the  out¬ 
put  tape  for  the  penetration ,  radial  penetration ,  and  energy  distributions 
of  the  knock-on  atoms  . 

BOX  10 

The  total  number  of  atoms  that  exited  through  a  lattice  face  is 
written  on  the  output  tape.  The  number  of  the  atom,  its  co-ordinates, 
its  velocity  components,  energy,  LB  (in  octal) ,  its  penetration,  radial 
penetration,  and  time  step  of  exit  are  written  on  the  output  tape  for 
each  atom  exiting  through  a  face . 

The  entire  box  is  then  repeated  for  each  of  the  six  faces  of  the 
lattice . 

BOX  11 

If  Jump  Switch  1  is  set,  a  jump  is  made  to  the  end  of  the  program, 
otherwise  a  general  lattice  dump  is  written  on  the  output  tape.  The  dump 
consists  of  the  atom  number,  co-ordinates,  velocity  components,  energy, 
and  the  LB  contents  (in  octal)  for  each  atom. 

SUBROUTINE  ROUND  (L) 

BOX  12 

This  subroutine  computes  the  site  closest  to  atom  L.  The  atom's 
co-ordinates  are  rounded  off,  and  the  fixed  point  co-ordinates  of  the 
atom  are  stored  in  the  IB  array  if  the  co-ordinates  represent  a  site. 
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If  the  fixed  point  co-ordinates  do  not  represent  a  site ,  the  six  nearest 
neighbor  sites  are  found.  The  site  that  is  closest  space-wise  to  the 
atom  is  chosen  as  the  site  for  atom  L.  Once  a  site  is  chosen,  a  RETURN 
to  the  main  program  is  executed . 

SUBROUTINE  INITIAL  (L) 

BOX  13 

If  atom  L  is  one  of  the  original  lattice  atoms ,  this  subroutine 
computes  its  initial  fixed  point  co-ordinates  and  stores  them  in  the  IB 
array.  If  atom  L  was  shot  into  the  lattice,  the  point  of  entry  is  stored 
in  the  IB  array .  A  RETURN  to  the  main  program  is  then  executed . 

SUBROUTINE  ECOUNT  (L) 

BOX  14 

The  energy  of  atom  L  is  inspected.  This  energy  will  fall  between 
two  of  the  values  on  the  KET  list,  say  the  sixth  and  the  seventh  values 
on  the  list.  Atom  L's  energy  is  then  greater  than  the  sixth  (ith)  value 
on  the  list.  The  histograph  counter  IC  (I  +  50)  is  then  increased  by  one. 
A  RETURN  to  the  main  program  is  then  executed. 

SUBROUTINE  DEPTH  (L) 

BOX  15 

The  penetration  (DEPTH)  is  computed  as  the  final  X  co-ordinate  of 
atom  L  minus  the  initial  X  co-ordimte.  A  value  of  I  is  determined  so 
that  DEPTH  lies  between  A  (I  -7.5)  and  A  (I  -6.5) .  Histograph  counter 
IC  (I)  is  then  increased  by  one . 
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The  radial  penetration  (PERDIST)  is  then  calculated  from  L's  initial 
and  final  positions.  A  value  of  I  is  determined  so  that  PERDIST  lies  be¬ 
tween  A  (I  -1)  and  A  (I) .  Histograph  counter  IC  (I  +  30)  is  then  increased 
by  one.  A  RETURN  to  the  main  program  is  then  executed. 

SUBROUTINE  HISTO  (NGROUPS) 

BOX  16 

The  largest  number  and  the  total  of  NUMBERS  (1  -NGROUpS)  are 
determined .  The  total  is  written  on  the  output  tape .  The  percentage  of 
the  total  is  written  on  the  output  tape  for  each  entry  on  the  NUMBERS 
list,  followed  by  a  number  of  X's  corresponding  to  the  percentage  this 
entry  is  of  the  largest  number  in  the  table,  i.e.  the  largest  number  will 
write  out  100  X's,  and  a  number  half  that  large  will  write  out  50  X's. 

A  RETURN  to  the  main  program  is  then  executed. 

PROGRAM  RON 

This  program  writes  on  the  output  tape  the  successive  positions  of 
every  atom  in  the  lattice  (or  that  was  shot  into  the  lattice) ,  if  the  suc¬ 
cessive  position  of  an  atom  is  significantly  different  from  the  previous 
position.  The  number  of  the  time  step  is  also  written  when  an  atom's 
co-ordinates  are  written.  After  writing  all  the  successive  positions  of 
an  atom,  the  program  commences  writing  the  positions  of  the  next  atom 
that  has  moved  significantly  until  all  atoms  have  been  considered. 
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